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An  exact  formulation  was  developed  for  the  prediction  of  the  internal  el^tnc  field  of  an  assembly  of 
Rayleigh  particles  exposed  to  an  incident  radiation  field  which:  (i)  satisfies  the  eneipr  conservation 
(ii)^  amounts  for  self  contribution  and  multiple  scattering  effectsini)  is  invanant  to  c^rfinate 
transformation  and  (iv)  is  particularly  useful  for  soot  diagnostics.  Based  on  this  formulafaon  a 
computer  code  (AGGL)  was  written  for  the  computation  of  exti^ion,  total  scattering  eftteiency 
factors,  and  differential  scattering  intensities  for  agglomerates.  The  computations  were  free  of 
error  and  resulted  in  tenfold  reduction  in  computational  time. 

Agglomerates  with  straight  chain  configuration  were  generated  in  edition 

Se  seeded  with  iron  pentacarbonyl  v^r.  Both  the  pnmary  particle  size  and  «>®  ^ 

were  found  to  be  strong  functions  of  theXvapor  concentration  a^ag^gate 

diffraction  measurements  revealed  that  the  particles  cxinsisted  predominantly  of  iron  oxide(Fe203). 

In-  situ  measurements  of  the  particle  dynamics  using  polariz^  imd  depolarized  light  ncatUm^ 
showed  that  the  signal  to  noise  ratio  was  significantly  improved  when  m^ure^nts  were  came 
out  in  the  cross  correlation  mode.  It  was  also  found  that  monitoring  of  the  particle  dynamics  was 
possible  with  sampling  times  as  low  as  26  nanoseconds. 
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1.  Summary  of  the  Work  Carried  out  During  the  Grant  Period 

The  work  carried  out  during  the  grant  period  may  be  summarized  as  follows: 

a)  An  exact  formulation  has  been  developed  for  the  prediction  of  the  internal  electric  field  of  an 
assembly  of  particles  exposed  tc^  an  incident  radiation  field  which:  (i)  satisfies  the  energy 
conservation  requirement  (ii)  accounts  for  self  contribution  and  multiple  scattering  field  effects  (iii) 
is  invariant  to  coordinate  transformation  and  (iv)  is  particularly  useful  for  soot  and  other 
particulate  diagnostics  at  the  ultraviolet  (UV)  part  of  the  spectrum.  This  part  of  the  work  has 
been  published  in  the  Journal  of  Physics  D:  Applied  Physics,  Vol.  27 , 2258-2270  (1994)  under 
the  title  “On  the  electromagnetic  scattering  and  absorption  of  agglomerated  small  spherical 
particles”. 

b)  Based  on  the  developed  exact  formulation  of  light  scattering  by  agglomerated  small  particles, 
partial  derivatives  of  scattering  quantities  were  derived  with  respect  to  refi’active  index  and  size 
parameter.  Computations  of  the  derivatives  could  be  conveniently  incorporated  into  scattering 
computer  code.  This  provided  for  an  efficient  computation  of  the  various  derivatives  needed  for  data 
inversion  purposes  as  opposed  to  the  usual  finite  difference  approximation  which  requires 
judicious  choice  of  step  sizes  and  additional  computational  time  and,  as  a  result,  is  subject  to 
inaccuracy  as  well  as  inefficiency. 

Specific  considerations  regarding  data  inversion  were  provided.  Because  experimental  errors  are 
inevitable,  it  is  suggested  that  an  error  response  or  inversion  uncertainty 

analysis  could  be  included  as  an  inherent  part  of  the  data  inversion  task.  A  complete  data 
inversion  procedure  should  consist  of  an  evaluation  of  derivatives,  inversion  of  measured  data  and 
estimation  of  imcertainty  resulting  finm  measurement  of  errors. 

The  pairing  problem  in  the  sense  of  deciding  which  quantities  to  measure  in  order  to  obtain  the 
most  accurate  and  reliable  results  was  also  posed  and  was  found  to  be  a  critical  and  important 
step. 

With  respect  to  the  inversion  scheme,  it  was  found  that  the  same  condusions^restrictions  apply  to 
agglomerate  geometries  other  than  the  straight  chain  one,  such  as  dusters(closely  packed  structure) 
and  randomly  branched  chain. 
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A  simple  inversion  algorithm  for  the  complex  refractive  index  was  developed  which  exhibits 
excellent  numerical  behavior  and  provides  a  prototype  for  developing  more  sophisticated  inversion 
algorithms.  Because  aU  further  data  inversion  of  light  scattering  data  inevitably  involves  inversion 
for  the  complex  refractive  index,  the  approach  developed  during  the  course  of  this  study  constitutes 
the  building  blocks  for  future  studies  fliming  at  complete  characterization  of  agglomerated  particle 
systems.  The  results  of  this  part  of  the  work  have  been  published  in  the  Journal  of  Physics  D: 
Applied  Physics  Vol.  28  ,  pp  2585-  2594(1995)  imder  the  title  “On  the  inverse  acattering  problem 
for  characterization  of  agglomerated  particulates”  partial  derivative  formulation” 

c)  Based  on  the  the  exact  light  scattering  formulation  developed  and  summarized  under  (a)  above 

as  new  computer  code  under  the  acronym  (AGGL)  was  written  for  the  computations  of  extinction, 
absorption,  total  scattering  efficiency  factors,  and  differential  scattering  intensities  for  agglomerates 
of  arbitrary  configurations.  Derivatives  with  respect  to  the  refractive  index  and  size  parameter  were 
also  internally  generated.  The  computations  showed  that  the  code  was  of  error 

and  resulted  in  tenfold  reduction  in  computational  time.  Details  about  the  developed  code  and  a 
listing  of  the  code  itself  is  presented  in  Section  VII  of  this  Report. 

d)  Since  one  of  the  objectives  of  the  study  was  to  develop  a  technique  for  the  determination  of  the 
dynamics  of  particulates  in  reacting  flows  a  system  was  developed  by  which  agglomerated 
particulates  of  specific  configurations,  namely  chain  hke  aggregates  could  be  produced.  The  system 
can  accommodate  both  mechanical  and  thermophoretic  particle  sampling  and  optical 
measurements  to  be  performed  as  functions  of  axial  and  radial  position  within  a  difiusion  flame 
system.  The  reactor  consists  of  a  diffusion  type  burner  constructed  of  stainless  steel.  The  fuel  used 
was  carbon  monoxide  seeded  with  iron  pentacarbonyl  vapor  (Fe(CO)5)  and  the  oxidizer  was 
ambient  air.  The  details  of  the  experimental  system  are  discussed  in  the  Section  IV  of  this  Report. 

e)  Using  the  experimental  system  developed  (see  (d)  above)  nanosized  chain— like  aggregates  were 
synthesized  in  an  iron  pentacarbonyl— carbon  monoxide— air  difiusion  flame  system.  The  size  and 
shape  of  the  chainlike  aggregates  and  primary  particles  were  characterized  in  terms  of  measurable 
moiphological  parameters  such  as  diameter  of  primary  particles  and  aspect  ratio.  The  influence  of 
operating  conditions  on  the  parameters  have  been  investigated  systematically  by  employing 
thermophoretic  and  transmission  electromicroscopy  (TEM).  The  morphological  parameters  were 
found  to  be  strong  functions  of  operating  conditions,  residence  time,  and  locations  within  flie  flame 
and  could  be  controlled  by  adjusting  the  concentrations  of  Fe(CO)5  seeded  in  the  flame  through 
the  variation  of  the  carrier  gas  flowrate  and  temperature  of  metal  additive. 

X— ray  diffraction  measurements  coupled  with  thermodynamic  equilibrium  composition 
predictions  were  used  to  identify  the  chemical  state  of  the  particulate  components.  The  results 
indicated  that  the  aggregates  formed  in  this  flame  consisted  predominantly  of  primary  particles  of 
Fe203.  In  addition  the  ph5^ical  and  chemical  mechanisms  of  formation  of  chain-like  aggregates 
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were  reviewed.  Our  observations  pointed  to  the  possibility  that  while  magnetic  forces  are 
necessary  for  providing  a  self  -  alignment  mechanism,  the  flame  temperature  may  also  be  critical 
the  formation  of  chain  like  aggregates.  The  results  of  this  part  of  the  work  were  presented  in 
Twenty  Sixth  Symposium  (International )  on  Combustion  and  Published  in  the  Proceedings  of  the 
Combustion  Institute  xmder  the  title  “Controlled  Combustion  Synthesis  of  Nanosized  Iron  Oxide 
Aggregates  “  and  co-authored  by  Z.  Zhang  and  T.T.  Charalampopoulos. 

f)  As  planned  in  the  original  proposal  an  optical  system  had  to  be  developed  that  would  allow 
monitoring  of  the  depolarized  correlation  function  corresponding  to  the  chain-hke  aggregates 
formed  under  flame  conditions.  The  optical  system  was  also  to  provide  for  monitoring  of  the 
polarized  correlation  function  from  the  same  set  of  aggregates.  This  particular  system  has  been 
developed  and  has  been  tested.  The  resultsobtained  so  far  have  shown  that  data  can  be  collected 
at  sampling  rates  as  low  as  25  nanoseconds.  Furthermore  it  has  been  shown  that  significant 
reduction  in  noise  levels  can  be  accomplished  when  the  measurements  are  carried  out  in  the  cross 
correlation  mode.  This  is  because  noisy  efiects  such  as  those  resulting  from  aflerpulsing  of  the 
photomultiplier  tubes  become  virtually  xmcorrelated  under  cross  correlation  conditions.  The  results 
of  this  part  of  the  work  are  presented  in  Section  V  of  this  Report. 

g) In  the  course  of  the  development  of  the  method  for  obtaining  the  dynamics/morphology  of  the 
aggregates  using  depolarized  light  scattering  the  use  of  heterodyne  dynamic  light  scattering  to 
infer  velocities  and  Brownian  diflPusion  coefficients  of  agglomerates  in  reacting  flows  was  assessed 
analytically.  The  analysis  related  the  variance  of  the  inferred  velocity  and  difiusion  parameters  to 
the  variance  of  the  normalized  correlation  function  of  three  detected  light’s  intensity  (or  photon  coimt 
)  fluctuations-  An  experimental  criterion  was  established  for  insiiring  the  accuracy  of  the  inferred 
results.  This  criterion  defined  the  range  of  angles  for  which  the  statistical  uncertainties  of  the 
inferred  results,  when  normalized  with  respect  to  the  inferred  parameters,  magmtudes  are  less 
than  or  equal  to  the  statistical  imcertainty  of  the  experimentally  determined  correlation  function. 
Under  suitable  assumptions,  this  analysis  is  extended  to  non-spherical,  cylndrically-- symmetric 
particles  (or  agglomerates) 
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IV.  Development  of  a  system  for  controlled 
combustion  synthesis  of  chain-like  aggregates 

1.  Introduction 


.  Knowledge  of  the  shape  and  size  of  particles  and  their  aggregates  formed  in  the  combustion 
system  is  significant  in  many  areas  of  research  and  practical  applications.  Information  about  the 
size 

and  shape  of  anisotropic  particulates  is  important  for  (i)  predicting  growth  and  oxidation  of 
particulates  in  combustion  systems,  (ii)  studying  the  agglomeration  mechanism  of  particulates,  (iii) 
the  efficient  control  of  fibrous  aerosols  in  the  industrial  hygiene,  (iv)  the  study  of  real  time  quality 
control  of  materials  synthesized  in  flame  reactors,  and  (v)  air  pollution  control  studies[l] 

The  techniques  used  for  determining  the  size  and  shape  of  particulates  in  flame  system  can  be 
divided  into  two  categories:  ex-situ  and  in-situ.  Both  techniques  possess  advantages  and 
disadvantages  and  are  complementary.  The  present  study  constitutes  the  first  step  towards  the 
development  of  an  in-situ  optical  method  by  which  the  optical  and  morphological  properties  of 
s^glomerated  structures,  as  well  as  their  dynamics  in  reacting  systems,  may  be  determined.  The 
potential  of  such  an  approach  and  its  implications  has  been  recently  discussed  in 
Charalampopoulos  [1].  Thus  the  objectives  of  this  part  of  the  study  were  to  develop  a  suitable 
experimental  system  for  the  synthesis  of  anisotropic  aggregates  and  to  determine  the  most 
appropriate  type  of  extractive  sampling  technique  without  comproimsing  the  integrity  of  the 
samples.Since  straight  chains  constitute  the  simplest  structure  and  most  weU  defined  form  of 

aggregates,  they  have  been  chosen  for  the  purpose  of  this  study. 

Numerous  experimental  studies  have  been  performed  to  sjmthesize  subnucron  size  particulates  and 
to  determine  the  shape  and/or  the  morphological  features  of  aggregates  formed  in  combustion 
processes  [2-9],  employing  both  ex-situ  and  in-situ  methods.  The  studies  show  that  the  structure  of 
aggregates  in  fiflirift  strongly  depends  on  the  technology  of  synthesis,  the  type  of  combustion  system, 
fuel,  fuel  equivalence  ratio,  and  operating  conditions.  Some  of  the  studies  specifically  investigated 
the  synthesis  of  chain-like  aggregates.  Two  techniques  were  used  in  the  synthesis  of  the  chain-like 
j^gregates.  With  the  exploding  wire  method  [2,3]  predominantly  chains  were  produced,  consisting 
up  to  thousands  of  primary  particles,  while  the  disadvantages  of  the  method  are  that  the  process  is 
not  continuous  and  is  difficult  to  control  in  order  to  produce  uniformly  dispersed  chains.  The 
combustion  synthesis  via  the  gaseous  phase  was  used  for  the  production  of  chain-like  aggregates  as 
well.  The  metal  additives,  such  iron  pentacarbonyl,  are  introduced  into  flames  where  decomposition 
and  chemical  reaction  take  place  to  form  metal  oxide  vapors.  Because  of  the  nature  of  the  chemical 
structures,  the  molecules  of  some  metal  oxides  ,  such  as  Fe203  [  5  ] ,  tend  to  form  chain  structures. 
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Then  the  coagulation  of  primaiy  particles  from  subsequent  nucleation  and  condensation  leads  to 
chain-like  aggregates. 

A  number  of  experiments  [  5,7  ]  have  been  conducted  to  synthesize  the  chain-like  aggregates  based 
on  the  above  combustion  technology.  The  chain  -like  aggregates  of  iron  oxide  particulates  were 
generated  by  seeding  a  C0/02/N2  flame  with  iron  pentacarbonyl  (  Fe(CO)5  )  vapor.  The  Fe{CO)5 
vapor  decomposed  at  above  300  C  to  form  iron  oxide  in  the  presence  of  oxygen.  Under  appropriate 
flame  conditions,  the  iron  oxide  particles  formed  chain-like  a^egates.  The  investigators  also 
characterized  the  morphology  of  chain-like  aggregates  employing  ex-situ  technique.  The 
characteristics  of  the  chain  aggregates  were  determined  for  the  samples  taken  in  the  post  flame 
region  for  different  concentrations  of  the  seeding  vapor  and  varying  02:C0  mixing  ratio.  The  studies 
found  that  both  diameter  of  primary  particles  and  length  of  chain-like  aggregates  were  controllable 
and  dependent  of  additive  concentration  [6,7]. 
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2.  Development  of  the  system 


In  order  to  fulfill  the  objectives  of  the  study,  the  experimental  sj^m  should  exhibit  tiie  following 
properties:  control,  reproducibility,  and  stability.  In  accordance  with  these  design  criteria,  the 
e^)erimental  system,  developed  fiir  the  controlled  combustion  synthesis  of  diam-like  aggregates  and 
fi>r  the  performance  of  the  ex-situ  characterization,  consists  of  a  supply  and  feeding  ss^tem,  which 
supplies  the  combustion  gases  and  seeding  vapor  to  the  burner  and  r^ulates  the  flowrates  of  all 
gases  in  the  combustion  process;  a  burner,  which  produces  the  flame;  sampling  system,  which 
extracts  the  particulates  fi:om  the  flame  to  facihtate  the  TEM  analysis;  and  driving  mechamsm 
systems,  which  controls  the  location  of  the  burner  with  respect  to  the  sampling  probe  and  optical 
apparatus.  A  schematic  of  flie  experimental  set-up  is  shown  in  F^ure  1. 

A  supply-feeding  system,  as  shown  in  Figure  2,  was  designed  to  deliver,  control  and  monitor  the 
gases  to  the  burner.  The  amount  of  each  component  fed  into  flie  burner  was  accurately  controlled 
and  monitored  by  this  system.  The  capabilities  of  the  system  include  the  ability  to  provide 
independent  variations  of  the  feeding  flow-rate  of  gases  through  a  network  of  compressed  gas 
cylinders,  flowmeters,  valves  and  tubing.  The  flowrates  of  critical  reactants,  CO  and  02  were 
monitored  by  means  of  the  thermal-type  digital  mass  flow-meters  [Hastings,Model  HFM  200]  with 
fluctuation  less  than  1  percent  of  the  full  scale  flow  rate  of  5  1/min.  Because  of  the  high  corrosion 
effect  of  iron  pentacarbonyl  vapor,  a  purging  system  is  incorporated  with  the  gas  supply  system  that 
the  whole  gas  supply  system,  including  the  burner  and  digital  flow  meters,  is  puiged  wnth  nitrogen 
after  every  run  of  the  flame.  The  purging  nitrogen  acts  to  remove  some  of  the  deposits  within  the 
reactant  supply  line  and  burner,  which  maintains  the  cleanliness  of  the  supply  and  feeding  system 
for  continuous  accuracy  and  provides  the  necessary  condition  of  reproducibility  and  control  of  the 
flame  experiments. 

Figure  3  illustrates  the  additive  feeding  system.  The  iron  pentacarbonyl  vapor  is  introduced  into 
the  CO-air  flame  by  diverting  a  small  ft’action  of  bulk  fuel  flow  through  a  column  of  liquid  Fe{CO)5 
stored  in  a  heated  metallic  stainless  steel  q^der.  By  monitoring  the  temperature  of  the  cylinder 
and  carrier  gas  flowrate,  partially  saturated  Fe(CO)5  vapor  is  transported  by  the  carrier  gas  and 
injected  into  the  flamft  throu^  the  inner  tube.  A  thermal-type  digital  mass  flow-meter  [  Hastings, 
Model  HFM200L  ]  is  used  to  monitor  the  flow  rate  of  carriage  gas  with  fluctuation  within  1  percent 
of  the  full  scale  flow  rate  of  300  cc/min.  The  concentration  of  iron  pentacarbonyl  vapor  introduced 
into  the  flame  was  regulated  by  controlling  the  temperature  of  the  heated  evaporate  cylinder  and 
the  carrier  gas  flowrate.  The  actual  mass  flowrate  of  iron  pentacarbonyl  supplied  for  a  given  carrier 
flowrate  and  cylinder  temperature  is  calibrated  by  measuring  the  condensed  mass  of  iron 
pentacarbonyl  vapor  during  a  measured  time  interval  into  a  condensation  coil  maintained  at  -15  C 
as  the  evaporator  outlet  gas  passes  through  the  coil.  Thus  the  concentration  of  iron  pentacarbonyl 
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vapor  delivered  to  the  bulk  fuel  flow  for  vaiying  carrier  flow  rate  at  a  given  cylinder  temperature  is 
known  prior  to  entering  into  the  flame.  The  calibrated  mass  flowrate  curve  along  with  calculated 
saturation  mass  flowrate  is  presented  in  Figure  4. 

In  order  to  minimize  the  condensation  and  deposition  of  the  iron  pentacarbonyl  vapor  on  the  wall 
along  the  inlet  line  during  the  flame  operation,  the  fuel  line  from  outlet  of  the  evaporate  cylinder  to 
the  inlet  of  the  burner  is  heated  to  an  av^age  temperature  92  C  and  thermally  insulated  from  the 
surroimdings.  Four  thermal  cx)uples  were  mounted  to  monitor  the  wall  temperature  at  various 
lcx;ations  along  the  line  fuel  supply  hne. 

The  burner  (  Figure  5)  used  in  this  study  is  a  concentric  stainless  steel  (316)  tube  burner  with  a 
diameter  of  the  inner  tube  1/4  inch  and  the  annular  interval  between  the  tubes  1/16  incii.  The  fuel, 
direcrted  to  the  surface  of  the  burner  through  the  central  tube,  mixes  with  throu^  diSusion  with 
nTidi7ftd^  whicih  is  either  obtained  frnm  ambient  room  air  or  directed  throu^  an  annular  r^on 
surrounding  the  fuel  tube  to  produce  a  nearly  conical  diffusion  flame.  A  ceramic  honeycomb  with 
fixed  above  the  burner  surface  was  utilized  to  straighten  the  gas  flow  streamline  and  to 

stabilize  the  flame. 

Two  different  sampling  systems  for  the  coUecrtion  of  particles  were  tested  in  the  course  of  this  study. 
One  is  a  mecdianic^al  suction  type  sampling  system,  the  other  one  is  a  thermophoretic  sampling 
system.  The  mechanical  extractive  type  sampling  system  is  depicted  in  Figure  6. 

’  Prior  to  burner  ignition  and  during  the  period  required  for  the  achievement  of  thermal  equihbrium 
between  the  burner  system  and  its  surroimdings,  a  small  amount  of  nitrogen  is  allowed  to  go 
thrniigh  the  probe  thrnngh  the  ‘purging’  line,  which  enables  the  sampling  system  to  be  kept  dean 
and  prevents  the  inlet  of  the  probe  frum  being  clogged  by  solid  particulates.  The  particle  sampling  is 
accomplished  by  establishing  a  vacuum  downstream  of  the  probe  and  consequently  drawing  the 
flflTTiA  gases  with  partides  into  the  probe.  The  by-pass  line  is  designed  for  the  selection  of  sampling 
rate  such  that  the  HiiniTmiTn  distuibance  of  the  flame  near  the  probe  inlet  is  observed,  which 
ensures  the  estabhshment  of  the  isokinetic  sampling  condition  [8].  When  the  desired  sampling 
condition  is  attained,  sampling  can  be  initiated  by  switching  the  two  three-way  valves  to  the 
‘sampling’  mode.  Upon  leaving  the  probe,  samples  pass  into  the  sampling  line  and  deposit  on  the 
carbon-coated  copper  grid  seated  on  a  Swagelok  in-line  filter  element. 

The  mechanical  suction  probes  designed  in  the  test  are:  a  water-cooled  probe,  and  a  uncooled  probe 
(as  shown  in  Figure  7).  The  original  overall  considerations  of  the  design  for  the  water  cooled 
sampling  probe  were  to  quench  frirther  agglomeration  and  chemical  reactions  within  tiie  probe  [  9  ]. 
It  consists  of  three  concentric  stainless  steel  tubes  (LD.of  0.125  inch,  O.D.  of  0.25  indi,  interval 
between  the  tubes  of  0.125  inch).  (Dooling  water  runs  between  the  two  annular  regions.  Samples  are 
withdrawn  throu^  the  center  tube.  A  thinner  probe  (I.D.  of  0.0785  inch,  O.D.  of  0.125  inch) 
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without  cooling  water  was  designed  to  minimize  disturbance  to  the  flame  while  chemical  reactions 
within  the  probe  were  of  less  concern. 

Since  continuous  agglomeration  within  the  probe  was  of  concern  in  this  study,  a  thennophoretic 
sampling  system  was  also  designed  and  constructed.  It  is  believed  [10]  tiiat  the  thennophoretic 
sampling  teduiique  is  well  suited  for  the  collection  of  particles  for  the  morphology  studies.  The 
principle  is  that  the  probe  provides  an  instant  cold  surface  in  the  hot  flame  to  establish  a 
temperature  gradient  which  drives  the  flame-bom  particles  toward  the  cold  surface  to  deposit  on  it. 
The  cold  surface  also  quenches  the  chemical  reaction  and  coagulations  of  particles  that  are  already 
captured.  Therefore,  the  thennophoretic  sampling  technique  eliminates  the  potential  agi^omeration 
of  particles[10,ll].  The  conventional  carbon  coated  copper  grid  [  Electron  Microscopy  Science,  Model 
CF200-Cu,  200mesh  ]is  used  as  a  cold  surface  to  catch  the  particulates.  The  carbon  coating  film  not 
only  is  stable  under  the  electron  beam,  comparing  with  other  coatings,  also  be  able  to  withstand 
relatively  high  temperature  as  well.  It  is  very  easily  oxidized  in  the  flame  environment,  however,  it 
can  be  exposed  in  the  flame  for  very  short  time  intervals.  A  slotted  plate  is  placed  between  the 
probe  and  flame  to  deflect  the  air  current  induced  by  the  motion  of  the  probe  and  to  block  the  heat 
transfer  from  the  flame  to  the  probe[  8  ] . 

The  thennophoretic  sampling  system  consists  of  a  double-acting  pnemnatic  cylinder  [BEMBA,  Model 
MRS-04-DXP]  of  3/4  inch  diameter  and  3  inch  stroke  that  is  driven  by  high  air  pressure  40  psi,  a 
solenoid  controlled  directional  4-way  valve,  and  a  computer  controlled  variable  time-delay  relay 
circuit  which  accurately  controls  the  position  and  motion  of  the  thennophoretic  probe  assembly.  The 
thennophoretic  probe  is  mounted  to  the  piston  shaft.  A  nitrid— based  permanent  magnet  is 
mounted  on  the  piston  to  operate  a  hall  effect  switch  (BIMBA  HSCX-04)  as  shown  in  Figure  8.  The 
BHE  switch  is  connected  in  series  with  a  D.C  power  supply  (9v)  and  a  oscilloscope.  During  the 
operation,  a  constant  current  is  passed  throu^  the  hall  sensor.  When  tiie  magnet  is  not  directly 
below  the  sensor,  the  current  distribution  will  be  uniform  and  no  potential  difference  exists  across 
tire  output.  When  the  magnet  is  directly  below  the  sensor,  it  disturbs  the  current  distribution, 
producing  a  potential  difference  across  the  output.  The  circuitry  amplifies  this  voltage  to  signal  the 
oscilloscope.  By  recording  the  time  variations  of  the  oscilloscope  signal  and  the  position  of  the  piston, 
the  position  and  movement  of  the  probe  can  be  precisely  determined.  A  piston  speed  of  0.86  m/s, 
and  a  transition  time  4+ 1ms  for  the  entry  of  the  probe  finm  the  edge  to  the  center  of  the  flame  are 
achieved  with  the  present  system.  The  computer  controlled  variable  time-delay  relay  circuit  also 
controls  the  sampling  time  which  varies  finm  75  to  200  ms,  depending  on  the  additive  concentration 
and  the  sampUng  position. 

The  precise  positioning  of  the  probe  with  respect  to  the  burner  surface  is  very  critical  in  order  to 
obtain  the  reproducible  measurements  of  size  distribution  varying  with  the  height  above  the  burner 
surface.  The  mechanical  components  of  the  thennophoretic  probe  control  system  and  mechanical 
suction  probe  system  are  mounted  on  a  rigid  rotational  optical  table.  To  provide  good  spatial 
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resolution  within  the  flame,  the  vertical  position  of  the  probe  is  fixed,  whereas  the  burner  moves  up 
flnd  down  through  a  burner  translation  mechanism.  A  pulse  generator  attached  to  the  gear  of  the 
burner  translation  mechanism  monitors  the  displacement  of  the  burner  such  that  the  hei^t  of  the 
probe  with  respect  to  the  burner  surface  can  be  accurately  controlled.  The  calibration  of  the 
translation  mechanism  determined  that  the  precision  for  measuring  the  actuate  burner  height  is 
within  1/10  of  a  millimeter.  The  burner  translation  mechanism,  as  shown  in  Figure  1,  is  similar  to 
the  one  used  by  Charalampopoulos  [  12  ].  The  precise  positioning  of  the  thermophotetic  probe  within 
the  flame  is  achieved  by  the  alignment  of  optical  arm  with  respect  to  the  laser  beam.  The  optical 
arm,  on  which  the  probe  is  mounted,  is  fixed  at  90  o  with  respect  to  the  incident  laser  beam  so  that 
the  probe  is  perpendicular  to  the  laser  beam  which  is  directed  and  focussed  at  the  central  axis  of  the 
burner.  The  carbon  grid  is  attached  on  the  tip  of  the  probe  and  its  position  is  zeroed  with  reference 
to  the  laser  beam.  The  precision  of  the  horizontal  positioning  of  the  probe/grid  assembly  is  within 
0.2  mm. 
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3*  Test  and  results 


The  diffusion  type  carbon  monoxide-air  flame  supported  on  the  system  was  tested.  The  tests 
showed  that  the  stable  seeded  CO-air  diffusion  flame  was  obtainable  for  CO  flowrate  ranging  from 
200  to  450  cc/min.  Flow  rates  lower  than  200  cc^nin  can  not  sustain  a  flame  and  flow  rates  higgler 
than  450  ccAnin  induce  flickering  in  the  flame.  At  a  total  fuel  flowrate  of  450  ccAnin  the  flame 
seeded  with  iron  pentacarbonyl  presented  a  most  stable  and  suitable  sampling  condition  for  the 
investigation  of  variations  of  the  particulate  morphology  with  the  operating  conditions,  with  the 
ceramic  honeycomb  fixed  at  40  mm  height  above  the  burner  siirface.  The  picture  of  this  flame  imder 
unseeded  and  seeded  conditions  is  shown  in  Figure  9.  The  flame  was  characterized  with 
temperature  measurement.  A  platinum  vs.  platmum-10\%  rhodium  thermocouple  with  a  bead 
diameter  of  0.3  mm  was  used  in  the  temperature  measurements.  The  position  and  alignment  of  the 
thermocouple  is  the  same  as  that  of  the  thermophoretic  probe.  Figure  10  shows  the  radiation 
corrected  radial  temperature  profiles  for  unseeded  CO/air  diffusion  flame  at  various  heists  above 
the  burner  surface.  The  maximum  fluctuation  for  the  temperature  measurements  occurred  at  35  mm 
height  and  was  found  to  be  2.7 \%  diffment  about  the  mean  temperature  value  of  three 
measurements. 

Experiments  for  synthesis  of  chain-hke  aggregates  and  sampling  particles  were  performed-  Chain¬ 
like  aggregates  were  produced  in  the  CO-air  diffusion  flame  seeded  with  Fe(CO)_{5}  vapor.  A  typical 
TEM  micrograph  is  shown  in  the  Figure  11  for  samples  taken  by  mechanical  suction  probes  and 
thermophoretic  probe.  The  image  of  the  iron  oxide  particles  appears  to  be  hexagonal  or  polyhedric 
whereas  the  image  of  the  ag^^ates  resembles  strai^t  chain-like  structure,  as  may  be  seen  from 
the  figure,  the  difference  of  the  results  was  significant  in  comparing  the  samples  collected  by  both 
mechanical  probes  with  thermophoretic  probe.  It  appeared  that  there  was  no  significant  difference 
in  the  primary  particle  diameter,  whereas  the  aggregate  morphologies  determined  by  these  two  ex- 
situ  techniques  were  highly  different.  The  size  and  shape  were  identical  for  the  samples  collected 
fium  the  water  cooled  and  the  non-cooled  probes.  Majority  of  the  aggregates  are  branched  chains 
and  consisted  of  hundred  of  primary  particles  per  aggregate,  yet  the  size  of  the  aggregates  could  not 
be  reduced  by  decreasing  the  concentration  of  Fe{CO)_{5}  vapor  and  did  not  vary  significantly  with 
location,  It  was  suspected  that  since  the  sizes  of  both  probes  are  comparable  with  that  of  the  flame, 
the  flame  might  have  been  severely  distorted  by  the  probe.  In  addition,  no  mechanism  exists  to 
prevent  additional  Agglomeration,  the  further  ag^omeration  might  have  occurred  within  the  probe. 
The  samples  collected  by  these  two  probes  probably  do  not  accurately  represent  the  true  morphology 
of  the  particulates  as  they  occur  within  the  flame.  However,  the  tiiermophoretic  sampling  technique, 
in  which  the  particulates  deposit  on  a  cold  surface  which  then  quenches  further  agglomeration, 
should  be  better  able  to  preserve  the  morphology  of  the  particulates  as  they  occur  within  the  flame. 
For  this  reason,  the  thermophoretic  sampling  system  was  utilized  for  the  determination  of  size  and 
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morphology  of  the  particles  and  their  aggregates  produced  in  the  experiments.  The  size  and  shape  of 
the  chain>like  a^^regates  and  the  primaiy  particles  were  characterized  in  terms  of  diameter  of 
primary  particles,  number  of  the  primary  particles  per  aggregate,  and  aspect  ratio  of  aggregates  of 
chain-like  aggregates  based  on  the  TEM  micrographs.  A  quantitative  data  on  the  primary  particle 
size  distribution  and  aggregate  size  distribution  as  a  function  of  location  in  the  jQame  were  obtained. 
The  variations  of  the  morphological  parameters  with  location  was  assessed  by  collecting  samples 
within  the  flame  along  the  centerline  at  the  heights  of  4,  15,  25,  35  mm  above  the  burner  surface. 
Figures  12  to  14  represent  the  frequency  distributions  of  the  diameters  of  primary  particles,  the 
number  of  primary  particles  per  aggregate,  and  the  aspect  ratio  of  the  strai^t  chains  as  determined 
by  aamplingr  and  TEM  photographs  at  the  operating  condition  of  total  fiiel  flowrate  450  aAnin  and 
carrier  flowrate  20  cc/min.  At  short  residence  times  (4mm  above  the  burner  surface)  the  straight 
chains  were  observed  with  ill  defined  shapes  that  quantitative  analysis  was  not  able  to  provided. 
As  the  height  increased,  the  particles  grow  and  the  time  available  for  aggregation  was  increased, 
resulting  in  longer  and  more  well  defined  chain-like  a^^'egates. 
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4.  Summary 

The  accomplishment  of  this  part  of  the  study  is  summarized  as  follows: 

1.  An  experimental  system  for  controlled  combustion  synthesis  of  particulates  has  been  developed  and 
tested  that  the  S5^tem  can  be  used  for  the  performance  of  both  the  ex-situ  morphol(^ 
measurements  and  the  in-situ  optical  measurements  of  the  determination  of  the  morphology  of  the 
anisotropic  particulates. 

2.Predominantly  chain-like  aggregates  were  synthesized  in  an  iron  pentacarbonyl  vapor  seeded 
co/air  diSusion  flame  supported  on  a  concentric  stainless  steel  tube  burner. 

3.  The  tests  showed  that  for  particle  morphology  studies,  thermophoretic  probe  system  may  be 
suitable  for  the  sampling  without  compromising  the  integrity  of  the  samples. 

4.  The  analysis  of  TEM  photographs  of  the  samples  collected  by  thermophoretic  sampling 
technique  provideda  qualitative  description  of  the  morphology  of  the  particulates  formed  in  the 
flame.  The  average  diameter  of  the  primary  particles  and  the  aspect  ratio  of  chain-like  aggregates 
increased  with  increase  of  the  residence  time. 
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Figure  1:  Experimental  setup 
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Figure  3:  Additive-feeding  system 
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Figure  4:  The  calibrated  aud  calculated  mass  flowrate  of  iron  pentacarbonyl  with  varying 
carrier  flowrate  at  heated  cylinder  temperature  2Z°C 
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Ceramic  honeycomb 
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Figure  5:  Concentric  tube  diffusion  burner 
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Figure  6:  Mechanical  suction  sampling  system 
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Figure  7:  Probes:  (a)  water-cooled  probe  (b)  uncooled  probe 
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Figure  8:  Pneumatic  actuator  and  thermophoretic  sampling  system 
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Figure  9:  CO-air  diffusion  flame  (a)unseeded  flame  (b)seeded  withFe(CO)s  ^^por 
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Figure  11:  The  TEM  micrographs:(a)  sample  taken  by  water-cooled  probe  (b)  sample  taken  by 
uncooled  probe  (c)sample  taken  by  thermophoretic  probe 
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Figure  12:  The  distributions  of  diameter  of  primary  particles  of  chain-like  aggregates  with 
varying  height  above  the  burner  surface 
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Aspect  Ratio 


Figure  13:  Tte  distribution  of  aspect  ratio  of  chain-like  aggregates  with  varying  height  above 
the  burner  surface 


y.  Use  of  Polarized  and  Depolarized  Dynamic  Light  Scattering  to  Determine  the 
Morphology  and  Dynamics  of  Flame-Borne  Agglomerates 

1.  Introduction 

The  objective  of  the  current  study  is  to  adapt  and  utilize  light  scattering  techniques  to  measure  both  the 
moiphological  (primary  particle  diameter,  aspect  ratio,  etc.)  and  dynamic  (Brownian  translational  and  rotational 
diffusion  coefficients)  properties  of  agglomerated  structures  occurring  within  flames.  Such  information  is  useful  for 
process  monitoring  of  materials  synthesis  using  flame  technology,  pollution  control  studies,  flame  emission 
calculations,  and  soot  formation  studies.  The  predominant  advantage  of  in  situ  measurements,  such  as  light 
scattering,  over  the  alternative  ex  situ  measurements,  such  as  the  extraction  of  particulate  samples  from  the  flame  for 
subsequent  examination  using  electron  microscopy,  is  that  the  in  situ  measurements  are  less  likely  to  disturb  the 
specimens  being  studied.  Additionally,  in  situ  measurement  techniques  are  generally  more  easily  automated. 

The  current  study  involves  performing  low  scattering  angle  (8  <  7®)  polarized  (W)  and  depolarized  (VH)  dynamic 
light  scattering  (DLS)  measurements  and  classical  depolarization  ratio  measurements  of  the  straight  chain 
agglomerates  of  iron  oxide  particles  occurring  along  the  central  axis  of  a  conical  iron-pentacarbonyl-seeded  C0-02 
diffusion  flame.  These  optical  measurements  are  supplemented  by  flame  temperature  measurements  using  a  type  S 
(Pt/Pt-10%  Rh)  thermocouple  and  an  ex  situ  determination  of  the  agglomerates’  morphology  distribution  obtained 
through  the  analysis  of  transmission  electron  micrographs  of  thermophorctically-obtained  samples. 

Utilizing  appropriate  theoretical  considerations  and  approximations  ([1],  [2],  [3]),  depolarized  dynamic  light 
scattering  data  can  be  combined  with  the  depolarization  ratio  and  either  polarized  DLS  data  from  the  same  scattering 
angle  or  depolarized  DLS  data  from  a  different  scattering  angle  to  infer  the  distributions  of  the  agglomerates’ 
Brownian  translational  and  rotational  diffusion  coefficients.  This  information  may  conceivably  be  combined  with 
the  measured  flame  temperature  and  appropriate  relationships  between  the  Brownian  diffusion  coefficients  and  the 
agglomerate  primaiy  particle  diameter  and  aspect  ratio  to  infer  the  mean  primaiy  particle  diameter  and  distribution  of 
aspect  ratios  of  the  agglomerates  at  the  measurement  location  within  the  flame.  Such  relationships  between  the 
Brownian  diffusion  coefficients  and  the  agglomerate  morphological  parameters  have  been  previously  theoretically 
obtained  for  particles  whose  characteristic  size  is  much  larger  than  the  mean  free  path  of  the  surrounding  gas  (the 
continuum  regime)  ([1],  [2],  [4]).  These  relations,  however,  are  not  valid  when  the  particle  size  is  of  the  same  cider 
(or  smaller)  as  the  gas  mean  free  path  [5],  a  condition  common  within  flames.  Because  the  non-continuum  effects 
on  these  diffusion  coefficient  relationships  are  not  accurately  known  at  the  present  time,  the  diffusion  coefficient 
information  obtained  from  the  optical  measurements  are  to  be  compared  with  the  ex  situ  morphology  information  to 
obtain  experimental  data  relating  the  Brownian  diffusion  coefficients  to  the  agglomerate  morphology.  Such 
information  would  contribute  to  the  development  of  a  non-continuum  model  or  empirical  correlation  relating  the 
diffusion  coefficients  to  the  morphological  parameters  so  that  the  morphological  informadon  can  be  directly  infmed 
from  future  DLS  experiments  without  the  need  for  ex  situ  sampling  and  analysis. 
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The  primary  motivation  for  the  use  of  dynamic  light  scattering  is  that  these  non-intrusive  measurements  exhibit 
little  or  no  depemlcnce  upon  the  typically  unknown  refractive  index  of  the  particles  being  measured  ([6],  [7])- 
Previous  DLS  flame  studies  (see  eg.  [6],  [8],  [9],  [10],  [11],  [12])  have  relied  exclusively  on  polarized  (W) 
measurements,  in  which  the  incident  and  detected  light  are  both  vertically  polarized  with  respect  to  the  (horizontal) 
scattering  plane,  which  is  defined  by  the  propagation  directions  (rf  tire  incident  and  detected  light.  These  polarized 
DLS  studies,  however,  typically  assume  that  the  scattering  particles  are  spherical  because  the  overwhelming 
dependence  of  polarized  DLS  data  on  the  isotropic  part  of  the  scatteters'  polarizability  tensor  permits  the  accurate 
determination  of  the  translational  diffusion  coefficient  only.  This  permits  the  extraction  of  an  equivalent-sjAere 
radius. 

This  investigation  extends  the  previous  DLS  flame  studies  by  incorporating  the  anisotropic  effects  obtained  from 
depolarized(VH)  dynamic  light  scattering  measurements.  Such  measurements  are  made  possible  only  by  the  recent 
improvements  in  the  minimum  sample  time  of  the  digital  correlators  (down  to  10-25  ns)  used  to  record  and  jnocess 
DLS  data.  The  rotational  decay  times  of  flame-bOTne  agglomerates,  which  can  range  from  10  to  100  ns  or  less,  are 
typically  too  fast  to  provide  a  meaningful  signal  when  the  scattered  light  is  processed  by  older  correlators,  whose 
minimum  sample  times  were  often  100  ns  or  longer.  With  the  additional  depolarized  information,  the  nonspherical 
nature  of  straight  chain  agglomerates  can  be  represented  through  the  determination  of  a  distribution  of  agglomerate 
aspect  ratios  in  addition  to  the  mean  primary  particle  diameter  based  upon  the  determination  of  both  the  Browrrian 
translational  and  rotational  diffusion  coefficients  from  the  DLS  data.  Since  straight-chain  agglomerates  are  the 
simplest  form  of  anisotropic  agglomerate  structure,  the  development  and  adaptation  of  light  scattering  methods  fw 
measuring  the  morphological  parameters  of  this  agglomerate  type  offers  a  step  towards  the  broader  objective  <rf 
developing  optical  techniques  to  measure  the  properties  of  arbitrary  structures. 
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2.  Apparatus 

The  apparatus  utilized  in  this  study  consists  of  a  concentric-tube  diffusion  burner,  which  generates  the  seeded  flame 
desoibed  above,  a  thermophoretic  sampling  system  for  extracting  ex  situ  samples  from  the  flame,  and  a  customized 
light  scattering  spectrometer  for  peifotming  the  optical  measurements. 

The  concentric-tube  diffusion  burner,  which  is  illustrated  in  figure  1,  is  constructed  of  three  concentric  stainless  steel 
tubes.  The  fuel  (CO)  flows  through  flie  central  1/2"  i.d.  tube,  emerging  to  combust  with  either  room  air  or  an 
annulus  of  O2  emerging  from  the  cylindrical  gap  surrounding  the  central  tube.  The  outer  tube  provides  an  outlet  for 
an  inert  shroud  gas,  such  as  N2,  which  can  be  used  prevent  entrainment  of  the  suirounding  air  into  the  flame.  The 
CO  fuel  is  seeded  with  iron  pentacarbonyl  (Fb(CO)5)  upstream  of  the  burner  by  diverting  a  portion  of  the  CO  flow 
through  a  temperature-regulated  cylinder  containing  200  ml  of  liquid  Fc{CO)5  und  subsequently  combining  the 
diverted  stream  with  the  main  CO  stream.  Bubbling  the  diverted  CO  flow  through  the  iron  pentacarbonyl  bath 
causes  partial  saturation  of  the  diverted  flow  with  the  volatile  Fe(CO)5.  The  tubing  between  the  Fe(CO)5  bath  and 
the  burner  is  electrically  heated  to  an  average  temperature  of  102°C  to  minimize  condensation  of  the  iron 
pentacaibonyl  on  the  tube  walls.  The  resulting  flame  is  stabilized  and  anchored  by  a  ceramic  honeycomb  (196  cells 
per  square  inch)  mounted  above  the  burner. 

A  thermophoretic  probe  ([13],  [14])  is  used  to  collect  samples  from  the  flame  for  an  ex  situ  determination  of  the 
distribution  of  agglomerate  morphological  parameters.  The  probe  actuator  inserts  a  room  temperature  20  x  4  x  0.4 
mm  stainless  steel  strip  into  the  flame  for  a  time  period  ranging  from  75  ms  to  200  ms.  The  temperature  gradient 
induced  in  the  flame  gases  between  the  room  temperature  probe  and  the  hot  flame  gases  propels  the  flame-borne 
particles  towards  the  probe  due  to  the  momentum  imbalance  of  the  molecules  impacting  the  particles  from  the  hot 
and  cold  sides.  This  phenomenon,  called  thermophoresis  ([15],  [16]),  causes  the  agglomerates  to  deposit  on  a  caibon 
coated  copper  grid,  which  is  attached  to  the  probe  strip.  The  particles  collected  on  this  grid  are  subsequently 
examined  using  a  transmission  electron  microscope  (TEM).  Disturbances  to  the  flame  arc  minimized  by  orienting 
the  probe  parallel  to  the  flame’s  streamlines  and  by  shielding  the  flame  from  air  currents  induced  by  the  ispdly 
moving  mechanism  by  placing  a  slotted  plate  between  the  actuator  mechanism  and  the  flame.  Theoretically,  this 
sampling  technique  does  not  bias  the  distribution  of  mori*ological  parameters  [17]  and  it  also  minimizes  the 
possibility  of  post-sampling  agglomeration  and  reaction.  It  is  thus  an  extremely  effective  available  ex  situ  sampling 
technique  for  morphological  studies  [14]. 

The  optical  system,  which  is  pictured  in  Figure  2,  is  a  light  scattering  spectrometer  of  fairly  standard  design  [18]. 
An  argon  ion  laser  beam  (from  a  Spectra  Physics  2085A-20  laser)  is  directedthrough  a  focusing  lens  LI,  an  aperture 
Dl,  and  a  polarizer  P  into  the  particle-laden  flame  to  provide  a  focused,  vertically-polarized,  cdierent  light  source  for 
scattering  by  the  agglomerates.  The  light  scatteredat  a  scattering  angle  0  from  the  direction  of  the  incident  beam  is 
again  polarized  by  the  analyzer  A  (either  in  the  vertical  plane  for  W  detection  or  the  horizontal  plane  for  VH 
detection)  before  being  focused  by  the  collector  lens  L2  onto  the  plane  of  aperture  D3.  A  photomultiplier  tube  PMT 
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(Thom  EMI  9863B/350),  whose  photocathode  is  immediately  behind  aperture  D3,  converts  the  photons  into  current 
pulses,  which  are  subsequently  discriminated  from  low  level  noise  by  a  preamplifier/discnminatOT  PAD  (Brocddiaven 
BIC  1(X)56).  The  PAD'S  ouq)ut  pulses  are  then  sent  to  a  Brookhaven  BI-9(X)0AT  digital  correlator  CORK,  which 
generates  the  following  correlation  function  [19]: 

Gixj)  =  ^  ^  ni  n'i-J  (/  =  1,2,3,...,M),  (1) 

i  =  1 

where  nj  and  nj-j  are  the  respective  numbers  of  photon  pulses  received  at  correlattH'  inputs  A  and  B  during  sample 
times  At  centered  at  times  tf  and  t;  -  Xj  (with  Xj  being  the  /th  delay  time),  N  is  the  number  of  prxxlucts  in  the 
sununation,  and  M  is  the  total  number  of  correlator  chatmels.  The  detection  system  ^rtures  D2  and  D3 
respectively  serve  to  define  the  effective  detector  area  and  the  width  of  the  scattering  volume  of  illuminated  particles 
seen  by  the  detector.  A  laser  line  filter  is  also  incorporated  within  the  detection  optics  to  eliminate  the  detection  of 
light  outside  the  wavelength  band  of  488  ±  1  nm.  The  high  voltage  power  supply  HVPS  shown  in  figure  2  is 
needed  by  the  PMT  so  that  the  photocurrent  pulses  can  be  amplified  to  detectable  levels  between  its  series  trf 
dynodes. 

In  addition  to  the  single-detector  autocorrelation  scheme  shown  in  figure  2,  a  two  detector  cross-correlation  technique 
is  also  considered.  As  shown  in  figure  3,  this  alternative  detection  system  uses  a  nonpolarizing  dielectric  cube  beam 
splitter  (Newport  05BC16NP.2)  to  direct  the  scattered  light  into  two  separate  photomultiplier  tubes.  The  resulting 
discriminated  signals  are  fed  into  the  A  and  B  inputs  of  the  correlator,  which  cross-correlates  them.  This  detection 
scheme,  as  desaibedby  Hiillies  [20],  minimizes  the  effects  of  afterpulsing  and  deadtime  on  the  correlation  function 
because  for  any  given  scattered  light  intensity,  the  specific  arrival  times  of  photons  detectedany  PMT  is  independent 
of  the  photon  arrival  times  at  the  other  PMT.  Because  both  afterpulsing  and  dead  time  effects  are  triggered  by  the 
arrival  of  individual  photons  at  a  detector  surface,  this  lack  of  correlation  between  the  photon  arrival  times  at  the 
different  detectors  greatly  reduces  the  afterpulsing  and  dead  time  effects  on  the  cross-correlation  function.  Some 
higher-order  afterpulsing  and  dead  time  effects  may  still  be  present,  however,  because  the  probability  of  detecting  a 
photon,  which  is  proportional  to  the  detected  intensity,  is  correlated  between  the  two  detectors.  The  main 
disadvantage  of  this  detection  method  is  that  tiie  scattered  light  intensity  detected  by  each  detector  is  less  than  half  cS 
the  intensity  detected  by  a  single  detector  in  the  autocorrelation  mode.  Thus,  the  experimental  duration  time  requited 
to  obtain  the  same  signal-to-noise  ratio  in  cross-correlation  experiments  as  obtained  in  autocorrelation  experiments 
must  be  increased  by  at  least  a  factor  of  four,  according  to  the  theoretical  uncertainty  results  of  Jakeman,  et  cd.  [21]. 
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3*  System  Tests  and  Preliminary  Results 

The  initial  tests  of  the  optical  system  involved  perfonning  measuiements  on  a  portion  of  the  incident  laser  light 
scattered  off  a  stationary  surface  into  the  detector  under  various  conditions.  These  static  light  experiments  were  used 
to  select  an  approfniate  supply  voltage  for  the  photomultiplier  tubes,  check  the  linearity  of  the  detection  system,  and 
detect  afterpulsing  and  dead  time  errors  in  the  measured  correlation  function. 

The  most  appropriate  supply  voltage  for  the  PMT’s  was  determined  by  measuring  the  detected  count  rate  as  a 
function  of  supply  voltage  at  a  constant  laser  intensity.  When  plotted  on  a  semi-log  graph,  as  in  figure  4,  the  count 
rate  versus  supply  voltage  curve  resembles  a  plateau,  with  the  count  rate  increasing  rapidly  with  increasing  voltage  at 
lower  voltages  and  leveling  to  a  neariy  constant  value  at  higher  voltages.  The  optimum  supply  voltage  is  chosen  as 
the  voltage  at  which  the  log  count  rate  versus  voltage  curve  flattens  [22].  This  insures  that  the  amplitudes  of  most 
of  the  photon-induced  current  pulses  from  the  photomultiplier  are  of  sufficient  amplitude  to  be  recognized  as 
photoelectric  pulses  by  the  preamplifier/discriminator  and  that  errors  and  fluctuations  in  the  supply  voltage  have  a 
minimal  effect  on  the  detected  photon  count  rate.  For  both  PMT’s  the  operating  voltage  was  determined  to  be 
approximately  1865  V.  Measuring  the  count  rate  at  various  laser  powers  subsequently  confirmed  an  jqiproximate 
linear  relationship  between  these  two  quantities,  in  accordance  with  theoretical  expectations  [22], 

In  order  to  check  optical  system  errors  in  generating  a  simple  correlation  function,  the  digital  correlator  was  used  to 
establish  the  time-based  correlation  function  of  the  detected  static  laser  light.  The  theoretical  correlation  function  of 
a  constant  intensity  source,  such  as  the  static  laser  light,  should  exhibit  no  variation  with  delay  time  (ie.  it  should  be 
a  flat  line).  Figure  5  depicts  the  normalized  autocorrelation  function  of  static  light  photon  counts  detected  by  a 
single  detector  with  the  optical  configuration  of  figure  2.  This  correlation  function  exhibits  systematic  noise  present 
from  0.025  ps  to  approximately  2  ps.  Initial  polarized  single-detector  DLS  measuiements  from  an  iron- 
pcntacarbonyl-seeded  CO/air  diffusion  flame  demonstrate  similar  noise  in  the  initial  channels  of  the  resulting 
correlation  function,  as  shown  in  figure  6.  A  comparison  of  figure  5  with  figure  6  reveals  that  the  noise  trends  are 
nearly  identical,  indicating  a  systematic  error  resulting  either  from  afterpulsing  and  deadtime  in  the  detection  system 
or  high-frequency  correlated  intensity  fluctuations  in  the  laser. 

In  order  to  determine  the  source  of  this  systematic  noise,  the  system  was  modified  to  incorporate  the  twodetector 
arrangement  of  figure  3.  The  resulting  normalized  static  light  cross-correlation  function  is  shown  in  figure  7.  This 
correlation  fimetion  was  measured  under  the  same  conditions  as  the  autocorrelation  function  of  figure  5  and  is  plotted 
on  the  same  scale  for  ease  of  comparison.  Figure  7  demonstrates  a  dramatic  reduction  in  the  systematic  noise  present 
in  the  0.025  ps  to  2  ps  range  in  the  static  light  correlation  function.  The  noise  is  also  noticeably  absent  from  the 
polarized  (W)  cross-correlation  function  depicted  in  figure  8,  which  was  generated  from  the  light  scattered  from  a 
seeded  flame.  This  result  implies  that  the  noise  in  the  single-detector  autocorrelation  function  was  created  by 
detection  system  effects,  such  as  afterpulsing  and  deadtime,  not  the  laser.  Because  the  0.025  ps  to  2  ps  delay  time 
range  is  crucial  for  accurate  determination  of  the  depolarized  DLS  correlation  function  from  the  flame-bome 
agglomerates  (due  to  its  expected  rapid  decay  time),  the  reduction  of  systematic  noise  within  this  range  outweighs  the 
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disadvantages  incuired  by  the  reduction  of  signal  resulting  from  splitting  the  scattered  beam  into  two  beams  prior  to 
detectioru 

Since  the  intensity  of  the  scattered  light^s  depolarized  component  is  typically  much  less  than  its  polarized  component 

(by  a  factor  of  10^^  or  less  for  agglomerates  occurring  in  the  presently-studied  flame),  the  conditions  necessary  for 
obtaining  a  depolarized  correlation  function  with  an  acceptable  signal-to-noise  ratio  are  much  more  extreme  than 
tiiose  required  for  polarized  DLS.  Based  upon  preliminary  measurements  in  conjunction  with  the  theoretical 
uncertainty  results  of  Jakeman,  et  al.  [21],  accurate  depolarized  DLS  measurements  require  an  incident  laser  power  on 
the  order  of  several  watts  with  experiment  duration  times  of  several  hours  (or  more).  Because  particle  deposition 
upon  the  burner  surface  and  stabilizing  hcmeycomb  is  observed  to  become  problematic  after  a  Fe(CX))5-seededC0-02 
flame  is  permitted  to  bum  for  longer  than  approximately  5  minutes,  an  accurate  depolarized  correlation  function 
cannot  be  obtained  by  collecting  the  data  continuously.  Therefore,  the  correlator  has  been  programed  to  accumulate 
correlation  function  and  count  rate  data  over  consecutive  30  second  intervals,  storing  the  data  accumulated  within 
each  interval  as  separate  datafiles.  This  permits  the  cessation  of  data  acquisition  every  5  minutes  to  permit  cleaning 
of  the  burner  surface  and  stabilizing  honeycomb.  The  resulting  multiple  correlation  functions  are  then  combined 
into  a  single  correlation  function  whose  effective  duration  time  is  the  sum  of  the  duration  times  of  the  individual 
correlation  functions.  This  technique  permits  the  effective  accumulation  of  DLS  data  for  arbitrarily  long  duration 
times  while  preventing  the  alteration  of  the  flame  conditions  due  to  excessive  particle  deposition.  This  technique 
also  facilitates  the  experimental  estimation  of  the  standard  deviation  of  each  point  on  the  combined  correlation 
function,  which  can  be  used  to  determine  an  effective  weighting  function  for  the  application  of  a  distribution  analysis 
routine,  such  as  Provencheris  CONTIN  ([23],  [24],  25]). 

Figures  8  and  9  depict  preliminary  polarized  and  depolarized  correlation  functions  of  photons  scattered  from  the  seeded 
flame.  The  total  CO  flowrate  was  0.50  1pm,  with  25  mlpm  diverted  through  a  room  temperature  (297  K)  bath  of 
iron  pentacarbonyl  prior  to  burning.  A  0.25  1pm  oxygen  flow  was  established  around  the  seeded  CO  stream  to 
facilitate  rapid  combustion  and  oxidation  of  the  iron  pentacarbonyl.  No  shroud  was  used  because  it  caused 
Instabilities  in  the  flame.  The  stabilizer  and  measurement  heights  were  45  mm  and  35  mm,  respectively.  The  light 
scattering  measurements  were  perfumed  at  a  scattering  angle  of  5®  with  an  incident  laser  power  of  approximately  1 
W  at  a  wavelength  of  488  nm.  The  total  duration  time  of  the  polarized  (W)  correlation  function  shown  in  figure  8 
was  5  minutes.  The  total  duration  time  of  the  depolarized (VH)  correlation  function  depictedin  figure  9  was  2  hours. 
The  polarized  DLS  measurement  resulted  in  a  reasonably  noise-free  correlation  ftmction  that  decays  to  its  baseline  at 
approximately  80  jjs.  The  noise  present  in  this  W  correlation  function  can  be  easily  reduced  by  extending  the 
duration  time.  The  VH  correlation  ftmction  exhibits  much  more  noise,  even  though  the  total  duration  time  was  24 
times  longer  than  the  duration  time  of  the  W  correlation  function.  Though  the  noise  level  is  too  great  for  a  reliable 
distribution  analysis  using  algorithms  like  CONTIN,  definite  trends  do  exist  in  the  depdarized  correlation  function: 
it  decays  from  an  initial  value  of  approximately  0.19  at  x  =  0.025  ps  to  its  baseline  at  approximately  x  =  0.5  ps. 
The  significant  temporal  variations  of  the  depolarized  correlation  ftmction  occur  entirely  within  the  range  of  delay 

34 


times  tiiat  would  have  been  significantly  plagued  by  afteipulsing  and  deadtime  had  the  single-detector  autocorrelation 
scheme  been  used 
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4.  Summary  and  Conclusion 

In  order  to  extend  currently  existing  techniques  of  measuring  the  morphological  parameters  of  flame-gencialed 
agglomerates  using  non-intrusive  light  scattering  techniques,  an  experiment  has  been  developed  to  measure  both  the 
polarized  and  depolarized  correlation  functions  of  laser  light  scattered  from  straight  chain  agglomerates  occurring 
within  an  mm-pentacarbonyl-seededCO-02  flame.  The  objective  is  to  obtain  distribution  functions  of  the  Brownian 
translational  and  rotational  diffusion  coefficients  from  these  correlation  functions  for  comparison  with  flame 
temperature  measurements  and  an  independent  ex  situ  determination  of  the  distribution  of  the  agglomerates' 
morphological  parameters.  Wth  this  information,  experimental  data  relating  the  diffusion  coefficients  to  the 
agglomerate  morphological  parameters  within  the  non-continuum  flame  environment  can  be  obtained.  Accurate 
knowledge  of  this  relationship  could  be  used  to  determine  agglomerate  morphology  distributions  directly  from 
dynamic  li^t  scattering  and  flame  temperature  measurements  without  the  need  for  intrusive  ex  situ  measurements. 
Preliminary  testing  of  the  apparatus  indicates  the  possibility  of  obtaining  both  optical  polarized  and  depdarLzed 
correlation  functions  from  a  flame  if  the  output  of  two  detectors  viewing  the  same  scattering  volume  are  cross- 
correlated  to  eliminate  unwanted  detection  system  errors,  such  as  afterpulsing  and  dead  time.  Under  the  conditions 
tested,  the  minimum  sample  time  of  the  correlator  used  to  process  the  data  must  be  shorter  than  100  ns  in  order  to 
resolve  the  temporal  variations  of  the  rapidly-decaying  depolarized  correlation  function.  Furthermore,  the  reduction 
of  noise  is  a  major  difficulty  in  measuring  the  depolarized  correlation  function  and  requires  high  laser  power  (several 
watts),  long  experimental  duration  times  (several  hours),  and  scattering  particles  exhibiting  strong  optical  anisotropy 

(depolarization  ratios  greater  than  10"^).  Continuing  work  is  underway  to  increase  the  signal-to-noise  ratio  of  the 
measured  depolarized  correlation  function  and  integrate  the  results  with  the  flame  temperature  and  ex  situ  morphology 
measuren^nts. 
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Figure  3 :  Schematic  diagram  of  two-detector  cross-correlation  detection  system. 
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Figure  6: 


Normalized  single-detector  autocorrelation  function  of  polarized  (W)  light 
scattered  from  an  iron-pentacarbonyl-seeded  CO/room-air  diffusion  flame. 
(CO  flow  rate  =  0.50  1pm;  carrier  CO  flow  rate  =  20  mlpm;  measurement 
height  =  30  mm  above  burner  surface;  stabilizer  position  =  42  mm  above 
iMiroer  surface) 
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Figure  8:  Normalized  two-detector  cross-correlation  function  of  polarized  (W)  light 
scattered  from  an  iron-pentacarbonyl-seeded  CO-O2  diffusion  flame.  (CO 
flow  rate  =  0.50  1pm;  carrier  CO  flow  rate  *  20  mlpm;  O2  flow  rate  = 
0.25  1pm;  measurement  height  =  35  mm  above  burner  surface;  stabilizer 
position  =  45  mm  above  burner  surface;  duration  time  =  5  minutes) 
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Normalized  Correlation  Function 


Figure  9:  Normalized  two-detector  cross-correlation  function  of  depolarized  (VH)  light 
scattered  from  an  iron-pentacarbonyl-seeded  CO-O2  diffusion  flame.  (CO 
flow  rate  =*  0.50  1pm;  carrier  CO  flow  rate  =  20  mlpm;  O2  flow  rate  = 
0.25  1pm;  measurement  height  =  35  mm  above  burner  surface;  stabilizer 
position  =  45  mm  above  burner  surface;  duration  time  =  2  hours) 
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VI.  Sensitivity  Considerations  for  the  Use  of  Heterod3Tie  Dynamic  Light 
Scattering  to  Infer  Velocities  and  Brownian  Diffusion  Coefficients  of 
Agglomerates  in  a  High  Temperature  Environment 

1.  Introduction  and  Motivation 

Experimeotal  detenninationof  the  morphology  of  flame-generated  particulates  within  a  flame  enviromnent  is  useful 
for  soot-formation  studies,  radiative  heat  transfer  computations  from  flames,  and  product  noonitoring  in  particulate 
synthesis  applications  [1].  Currently  the  available  in  situ  methods  for  particulate  morphology  diaracterization,  such 
as  the  fractal  analysis  technique  [l],are  limited  because  the  data  inversion  requires  knowledge  of  the  particle  refractive 
index,  which  in  most  cases  is  unknown.  Thus,  the  objective  of  the  ongoing  research  is  to  adapt  flie  technique  rf 
depdarized  dynamic  light  scattering  (DDLS),  which  currently  has  widespread  use  in  the  characterization  rf 
mactomolecules  in  solutions  [2,  3,  4]  and  is  independent  of  the  particle  refractive  index  for  particles  that  are  small  in 
comparison  to  the  incident  laser  wavelength,  to  the  measurement  of  the  dynamics  of  chain-like  agglcnnerates  within 
a  reacting  flame  environment.  Independent  extractive  sampling  and  subsequent  transmission  electron  microscope 
(TEM)  analysis  of  the  flame-generated  particulates  would  then  provide  valuable  data  on  the  relationship  between  the 
structure  and  non-continuum  rotational  dynamics  of  cylindrical  particulates.  Such  information  could  be  used  to 
confirm  any  proposed  theoretical  models  and  to  infer  particulate  structure  information  from  experimentally 
determined  rotational  dynamics  data. 

In  the  case  of  zero-angle  DDLS  measurements,  the  incident  laser  beam  is  sent  through  the  flame,  which  is  located 
between  crossed  polarizers,  and  straight  into  the  detector.  Ideally,  the  only  light  actually  reaching  the  detector  would 
be  the  depolarized  portion  of  the  scattered  light  caused  by  the  effective  optical  anisotropy  resulting  from  the  shape  of 
the  particulates.  The  autocorrelation  function  of  the  resulting  detected  intensity  fluctuations  would  then  contain 
rotational  diffusion  information,  with  no  dependence  on  either  the  translatimial  diffusion  coefficient  or  the 
interparticle  form  factor  [1,2, 4].  Because  the  scattering  volume  for  the  zero-angle  measurement  intersects  the  entire 
width  of  fire  flame,  the  effects  of  velocity  and  temperature  gradients  within  the  diffusion  flame  should  be  ctuisidered. 
If  these  gradients  significantly  affect  the  measured  correlation  function,  independent  measurements  of  radial  bulk 
particle  velocity,  as  well  as  temperature,  may  become  necessary  for  accurate  data  inversion.  Since  bulk  panicle 
velocity  can  theoretically  be  inferredfiom  the  results  of  a  heterodyne  DLS  experiment,  in  which  the  scattered  light  is 
mixed  with  an  unscattered  beam  (the  local  oscillator)  prior  to  detection  to  establish  a  mutual  beating  signal  [2,  3,  5, 
6],  the  remainderof  this  section  addressesthis  subject 
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2.  Experimental  Considerations 

Particles  that  are  suspended  within  a  fluid  (solvent)  medium  typically  undergo  a  variety  of  motions,  both 
translational  and  rotational,  that  generaUy  affect  the  temporal  variation  of  light  scattered  from  the  particle.  The 
particle  motion  is  often  governed  by  both  deterministic  and  random  forces.  Deterministic  particle  translation,  which 
is  here  characterized  by  a  bulk  or  average  velocity  V,  is  frequently  the  result  of  a  variety  of  phoretic  forces,  such  as 
tiiermophoresis,  photophoresis,  diflusiophoresis,  or  electrophoresis,  in  which  a  constant  net  force  is  ^{riied  to  the 
particle  and  causes  it  to  drift  in  the  force’s  direction.  Deterministic  translation  also  results  from  Stoke's  drag  induced 
on  the  particle  by  a  nonzero  solvent  velocity.  Additionally,  however,  the  dynamics  erf  microscopic  particles  arc  often 
significantly  affected  by  Brownian  motion,  which  is  caused  by  the  random  diermal  variations  of  the  momentum 
imparted  to  the  particle  by  the  solvent  molecules.  In  its  simplest  form,  Brownian  translation  is  desoibed  by  the 
following  diffusion  equation  [2]: 

^  Gs(R,t)=DV2Gs(R,t).  (D 

where  Gs(R,t)d^R  is  the  probability  of  finding  the  particle  within  the  neighborhood  (pR  of  position  R  at  time  t. 
The  proportionality  constant  D  is  called  the  Brownian  translational  diffusion  coefficient.  Rotational  Brownian 
motion  is  characterized  by  a  similar  probabilistic  diffusion  equation  containing  a  Brownian  rotational  diffusion 
coefficient  Dr  [2].  Autocorrelating  the  intensity  fluctuations  of  light  scattered  from  these  particles  at  different  times 
permits  the  experimental  determination  of  the  magnitudes  of  these  diffusion  coefficients,  which  can  then  be 
theoretically  related  to  the  particles’  size  and  aspect  ratio.  Additionally,  if  the  scattered  light  is  mixed  with  an 
unscattered  reference  beam,  a  beat  signal  is  established  between  the  two  beams  that  provides  information  about  the 
bulk  particle  velocity.  These  measurements  are  perfcMmed  using  the  technique  of  dynamic  light  scattering  (DLS). 

In  a  typical  dynamic  light  scattering  (DLS)  experiment,  polarized  laser  light  is  projected  through  a  suspension  of 
particulates.  The  resulting  scattered  light  then  passes  through  another  polarizer  (the  analyzer)  and  is  focused  onto  a 
photomultiplier  tube  (PMT).  The  PMT  output  is  then  processed  and  correlated  by  a  preamplifier/disciiminator  and  a 
digital  autocorrelator  to  provide  information  about  the  temporal  fluctuations  of  the  scattered  light,  from  which  the 
particle  dynamics  can  be  inferred.  Figure  1  depicts  the  scattering  geometry  used.  The  scattering  plane  is  the  plane 
defined  by  the  directions  of  the  incident  and  detected  scattered  beams,  whose  wave  vectors  are  ki  and  kf,  respectively. 
The  angle  between  ki  and  kf  is  the  scattering  angle  9.  The  scattering  vector  q  =  ki  -  kf  thus  has  the  magnitude 
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where  \  is  the  wavelengtib  of  the  incident  light  and  n  is  the  refractive  index  of  the  surrounding  medium.  The  basis 
vectors  characterizing  the  polarization  of  the  incident  and  scattered  fields  are  defined  to  be  parallel  (horizontal)  and 
perpendicular  (vertical)  to  fte  scattering  plane.  Typically,  the  incident  laser  light  is  vertically  polarized  while  the 
detectedlight  is  either  vertically  (for  W  detection)  or  horizontally  (for  VH  detection)  polarized.  The  experimental 
apparatus  is  discussed  in  further  detail  by  Ford  [7]. 

The  photomultiplier  tube  detects  the  individual  photons  of  the  scattered  light,  which  are  then  discriminated  from 
background  noise  using  a  preamplifier/discriminator.  These  photon  pulses  are  then  counted  and  processed  by  a 
digital  autocorrelator,  which  evaluates  the  following  correlation  function  [8]; 

<n(0)n(tj))  =  ^  ^  ^  ni  ni-j  ,  (j  =  1,  2, ...,  M)  (3) 

i  =  1 

where  U]  is  the  number  of  photons  detected  during  sample  time  At  centered  at  time  tj,  nj  is  the  number  of  photons 
detectedduring  At  centeredat  time  tj  -  Tj,  where  Xj  is  the  jth  delay  time,  and  M  is  the  number  of  correlator  channels. 
Neglecting  detection  system  effects  such  as  afterpulsing,  dead  time,  and  dark  current,  the  digital  autocorrelation 
function  is  related  to  the  normalized  second  order  electric  field  autocorrelation  function  (for  sample  times  much  less 
than  the  correlation  function's  characteristic  decay  time  [9])  by  the  following  equation  [3, 9, 10]: 

<n(0)n(x)>  =  <n(0))2|l  +  j^g(2)(x)  -  1]|  +<n(0))u(x),  (4) 

where  (n(0))  is  the  average  number  of  photons  counted  during  sample  time  At,  Nc  is  the  number  of  detected 

coherence  areas  [3, 6],  u(t)  is  a  function  whose  value  is  unity  for  t  =  0  and  zero  otherwise,  g®(T)  is  the  nOTtnalized 
second  order  electric  field  autocorrelation  function  with  the  definition 


{mO)|2|E(x)|2) 

(|E(0)i2)'2 

and  E(t)  is  the  electric  field  at  a  point  on  the  detector  at  time  t.  Equation  (4)  is  obtained  under  the  assumption  that 
the  first-oidertemporal  ctHrelation  function  of  the  energy  flux  within  each  coherence  area  is  equal  for  all  cdierence 
areas.  Additionally,  use  of  the  definition  of  g^2)(x)  in  equation  (4)  also  relies  on  the  validity  of  the  assumption  that 
the  electric  field  is  uniform  within  each  coherence  area.  The  first  term  on  the  right  side  of  equation  (4)  contains  the 
usable  signal  from  which  the  dynamics  of  the  scatterers  can  be  extracted.  The  second  term  is  a  shot  noise  term 
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resulting  from  the  Poisson  nature  of  the  photodetection  process.  When  the  cOTelation  function  is  defined  in  terms  of 
the  detected  photocurrent,  the  shot  noise  term  is  proportional  to  the  dirac  delta  function  S(x)  [11].  However,  since 
the  correlation  function  used  in  photon  counting  applications  is  formulated  in  terms  of  the  number  of  photons 
counted  during  a  finite  sample  time  At,  the  shot  noise  term  in  this  case  is  proportional  to  the  finite  step  function 
u(T)[9, 10]. 

For  a  heterodyne  experiment  where  the  intensity  of  the  local  oscillator  l]jj  is  much  larger  than  the  average  intensity 
is  of  the  scattered  field,  g^^%)  becomes  [3, 7]: 


2 

g(2)(x)  =  1+2  Re[Ii(T)  expCiiDLOt)],  (6) 

ilo 

where  Ii(x)  s  (Es*(0)Es(x))  is  the  first-order  scattered  electric  field  autocorrelation  function,  cdlo  is  ii*®  angular 

frequency  of  the  local  oscillator,  and  is  the  heterodyne  mixing  efficiency  [3].  For  a  monodisperse  suspension  of 
spherical,  optically  isotropic  particles,  Penner,  et  d.  [12]  has  evaluated  Ii(x),  under  the  assumptions  of  particle 
independence  and  a  spherically  symmetric  Gaussian  incident  intensity  distribution,  to  be 


Il(^)  =  is  exp(iq-VT)  expC-q^Dx)  exp  - 


vVl 

2(1)2 


exp(-i(DLi;), 


(7) 


where  V  is  the  bulk  velocity  of  the  particles,  D  is  the  particles'  Brownian  translational  diffusion  coefficient,  (U  is  the 
radius  of  the  scattering  volume,  and  (ul  is  the  angular  frequency  of  the  scattered  (and  incident)  light 
Substitution  of  equation  (7)  into  equations  (4)  and  (6)  gives  (for  (olq  =  (OL): 


<n(0)n(T))  =<n(0))2 


g2  ig  - 

1  +  2  ^  ~^cos(q-Vi:)  expC-q^Dr)  exp  -  — — 
Nc  iLO  ,  2(0^  , 


+  ^(0)>u(x).(8) 


When  (o/V  »l/q2D,  the  last  exponential  factor  is  negligible  because  its  decay  is  much  slower  dian  the  decay  of  the 

other  terms  in  equation  (8).  Assuming  that  both  the  baseline  (n(0))2  and  the  constant  in  front  of  the  cosine  factor 
are  experimentally  known  and  that  the  Gaussian  term  decays  much  slower  than  the  other  terms,  the  heterodyne 
correlation  function  can  be  normalized  into  the  form: 

g(x)  =  cos(At)  exp(-rt),  (forT  >  0)  (9) 
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where  A  =  q  V  andF  =  qpD. 


3.  Conditions  for  Particle  Velocity  and  Size  Inference  from  the  Heterodyne  Correlation 

Function 

An  analysis  of  the  variance  of  the  constants  A  andP  inferred  from  normalized  correlation  function  datagmC^iX  which 
is  obtained  atM  discrete  delay  times  Ti  (i  =  1,2,.  provides  the  basis  of  a  criterion  that  can  be  used  to  select  a 
range  of  scattering  angles  that  permits  the  inference  of  A  and  F  with  a  percent  uncertainty  equal  to  or  less  than  the 
uncertainty  of  the  normalized  taw  data  gm(^i).  The  data  inversion  is  assumed  to  take  the  form  of  a  least  squares 

analysis;  finding  the  parameters  A  andP  that  minimize  the  fimction  S?  s  EfemCi^i)  -  8(''i)l^-  following 

theoretical  analysis  utilizes  a  small  perturbation  approximation  to  determine  the  effects  of  the  statistical  uncertainty 
of  the  measured  correlation  function  gmCt)  on  the  statistical  uncertainty  of  the  least  squares  estimations  of  the 

parameters  A  and  T. 

To  facilitate  this  analysis,  gmC^j)  decomposed  as  follows: 


gm(^j)  =  g(^j)  +  o(^j)’ 

where  ^Tj)  is  the  value  of  the  theoretical  normalized  correlation  function  at  delay  time  Tj  and  o(Xj)  rs  the  deviation 
of  the  cwrelation  function’s  measured  value  from  its  theoretical  value  at  xj.  This  decomposition  assigns  all  of  the 
statistical  uncertainties  associated  with  the  measurement  to  the  random  function  o(Xj).  The  theoretical  correlation 
function  g(Xj)  is  therefore  a  deterministic  function. 

Utilization  of  a  least  squares  analysis  to  determine  the  unknown  parameters  A  andF  from  a  set  of  raw  data  involves 
detertnining  the  values  of  A  and  F  that  minimize  the  function 


,  f  2 

c2g  ^[gm(tj)  *  g(xj)] 
j  =  1 


[«m(tj)  -  cos(Axj)  cxp(-Fxj)] 


2 


j  =  1 


(11) 


where  Mis  the  number  of  correlator  channels.  In  this  case  the  values  of  A  andF  that  emrespondto  the  minimum  of 
c2  are  found  by  solving  the  equations  formed  by  setting  the  partial  derivatives  of  with  respect  to  A  and  F  equal  to 
zero: 
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gm(tj)  -  cos(ATj)exp(-rTj)]  xj  sin(Axj)  exjK-Fxj)  =  0 


(12a) 


j  =  1 


aid 


-  cos(ATj)exp(-rxj)]  Xj  cos(Axj)  exp(-rxj)  =  0 


(12b) 


j  =  1 


For  the  case  where  gm('tj)  follows  the  form  of  the  theoretical  curve  exactly  (gnif^j)  =  cos(AoXj)  exp(-roXj),  j  =  1, 
2, ....  M)  the  computed  least  squares  parameters  take  the  values  A  =  Aq  and  T  =  Fo-  When  gm('tj)  deviates  from  the 
theoretical  curve  by  the  random  quantity  of(Xj)  (ie.  gm(i^j)  =  cos(AoXj)  exp(-roXj)  +  o(Xj),  j  =  1,  2,  ....  M)  the 
computed  least  squares  parameters  diiferfrom  the  values  Ao  andF©  by  SA  and  S’ .  Thus,  for  this  more  general  case 
A  =  Ao+  SAandr  =  ro  +  S'.  By  substituting  gin(‘'j)  =  cos(AoXj)exp(-roXj)+ o(Xj),  A  =  Ao  + 6a,  andr  =  r, 

+  8r  into  equations  (12)  and  assuming  that  the  orderof  magnitudes  0[a(tj)]  =  0[5A]  =  OfS"  ] «  1 ,  equations  (12) 
approximated  by  their  first  m-derterms  in  o(xj),  8a,  and  S' ; 


o 

are 


o(xj)  +  Xj  [  S’  cos(AoXj)  +  8Asin(AoXj)]  exp(-roXj)}  xj  cos(AoXj)  exp(-roXj)  -  0 

j  =  1 


(13a) 


sod 


{  o(Xj)  +  Xj  [  S’  cos(AoXj)  +  8Asin(AoXj)]  exp(-roXj)}  xj  sin(AoXj)  exp(-roXj) 


=  0. 


j  =  1 


(13b) 


The  coefficients  of  the  perturbation  parameters  8A  and  ST  in  the  above  equations  are  linear  combinations  of  the  sums 


^  Xj2  cos(2AoXj)exp(-2roXj)  ,  Xj2  sin(2AoXj)  exp(-2roXj)  ,  and  2.  exp(-2roXj)  . 

j  =  1  j  =  1  j  =  1 

The  first  two  of  these  sums  are  the  real  and  imaginary  parts  of 
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(14) 


while  the  last  sum  is 


I  5  ^  Xj2  expI-2(ro  -  iAo)i:j]  , 

j  =  1 


Xj2  exp(-2roXj)  =  t 

j  =  1 


(15) 


Solving  equations  (13)  fw  5A  and  S'  and  expressing  the  results  in  terms  of  |  and  |o  gives  the  following 
relationships: 


[|o+  Re(D]  t  XjO(Tj)sin(AoXj)exp(-roXj)  -  Im(|)  xicO(Xk)cos(AoXk)exp(-r qX^) 

8a=2 - - 


k  =  1 


lt|2  -  tt? 


and 


(16a,b) 


[|o-Re(D]  Tjo(tj)cos(AoTj)exp(-roTj)  -  Ini(^  XkO(^k)sin(AoTk)e>^T otk) 


k  =  1 


By  computing  the  statistical  averages  of  the  squares  of  equations  (16),  first-order  pertuibation  approximations  of  the 
variances  of  the  least-squaies-inferred  quantities  A  and  T  are  obtained.  Squaring  the  expressions  for  8A  and  8" 
produces  terms  involving  products  of  the  form  otXj)o(xk);  therefore,  evaluation  of  the  variances  of  A  and  F  requires 
knowledge  of  the  covariance  matrix  Ck)v[ofTj),cj(Tk)l  =  {’cy(Tj)cr(Tk)X  If  mean  value  of  the  measured  correlation 
function  equal  to  g(xp  [=  cos(AoTj)  exp(-roXj)]  for  all  delay  times  Xj,  then  Cov[o(Xj),0(Xk)] 

=  Cov[gi„(Tj)^m(i^k)]- 

If  the  covariance  matrix  is  approximated  as  (a(Xj)o(Xk))  =  Cov[gixi(Xj),gin(^k)]  “  ^  where  is  the  vanance 
of  the  normalized  experimental  data,  which  is  assumed  to  be  independent  of  the  delay  time  Xj,  then  the  variances 
computed  using  the  expressions  for  SA  and  ST  given  by  equations  (16)  are 
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Var(A)  s  <(8A)2>  =  2  o2  ^  ^ 

lt|2  -  lo^ 


(17a) 


and 


Vai(r)s<(a')2)  =  2 


-,lo-Re(D 


(17b) 


Nondimensionalization  of  these  variances  produces  the  following  expressions: 


VaifA)  1  7o  +  Re(v) 
o2a2  2rt^ 


and 


VaifT)  2  Vo-ReCy) 
o^r2  (ix)2  vo^  -  *y'^ 


(18a,b) 


where 


(At)2 


Tj2  exp[-2x(r  -  i2u)Tjl 


=  1 


and 


t 


Tj2  exp[-23(TTj] 


j  =  1 


(19a) 


(19b) 


In  equations  (18)  and  (19)  the  dimensionless  parameters  r,  x,  andTj  are  defined  as 


r 


period 
decay  time 


2£o 

Ao  ’ 


_  sample  time 
period 


ApAt 
27t  ’ 


and  Ti. 

J  sample  time  at 


If  the  correlation  function  is  obtained  using  a  linear  correlator,  in  which  the  delay  time  is  an  integral  multiple  of  the 
sample  time,  ie.  tj  =  jAt  (orTj  =j),.then  equations  (19)  reduceto 
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(20a,b) 


s  s  exp[-2x(r  -  i2n)] 


yo=  ^  sJ  . 
j  =  1 


where 

and  So  s  exp(-2xr). 


The  series  expressed  in  equations  (20)  are  simply  the  second  derivatives  of  geometric  series.  Therefore,  both  series 
sum  to  the  following  closed-fcMm  expressions: 


V  =  s 


1  +  s  -  (M  +  1)^  sM  +  (2m2  +  2M  -  n  s^+l  -  8^+2 

(1  -s)3 


(21a) 


aid 


■yo  =  so 


1  +  So  -  (M  +  1)^  SqM  +  (2m2  +  2M  -  1)  SqM+I  -  Sq^+Z 

(1  -  So9 


(21b) 


which  reduce  to  the  following  as  the  number  of  conelator  channels  M  ^proaches  infinity: 

■»  =s(l.-Lsl  3^  •  (22a,b) 

(l-s)3  (l-so)5 

The  above  expressions,  which  relate  the  variances  of  the  least-squares-mferred  parameters  A  and  F  to  the  variance  of 
the  measured  normalized  correlation  function,  are  useful  for  determining  the  estimated  statistical  uncertainties  of  A 
and  r.  The  dimensionless  variance  ratios  on  the  left  sides  of  equations  (18)  are  the  ratios  of  the  squares  of  the 
inferred  parameters' percent  uncertainties  to  the  square  of  the  normalized  experimental  correlation  function’s  percent 
uncertainty  (as  T  — »  0).  These  variance  ratios  provide  guidelines  for  the  determination  of  appropriate  experimental 
parameters  to  produce  accurate  results.  All  subsequent  results  jnesented  in  this  paper  are  obtained  using  the 
expressions  for  y  and  yo  given  by  equations  (21)  and  (22),  which  are  valid  only  when  the  correlation  function's 
measured  delay  times  arc  uniformly  spaced.  Extension  of  these  results  to  multi-tau  correlation,  in  which  the 
measured  delay  times  are  exponentially  spaced,  can  be  accomplished  by  evaluating  the  sums  in  equations  (19)  with 
appropriately  spaced  Tj  values. 
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For  [Hesent  purposes,  a  maximum  "acceptable"  limit  of  uncertainty  is  defined  to  be  the  same  percent  uncertainty  as 
the  uncertainty  of  the  initial  normalized  data.  This  limit  is  set  in  order  to  insure  that  random  errors  in  the  raw  data 
do  not  produce  excessive  errors  in  the  inferred  results.  Thus,  setting  the  normalized  variance  ratios  expressed  above 
equal  to  unity  defines  the  boundaries  of  the  acceptable  region  of  operation  for  the  measurement  system.  These 
boundaries  are  shown  for  various  numbers  of  correlator  chaimels  M  in  figures  2  and  3.  All  points  in  the  plane  below 
the  curves  rejnesent  acceptable  operating  conditions  according  to  the  stated  criterion.  These  results  are  quite  general 
for  the  discrete  sampling  of  any  signal  with  the  form  of  equation  (9).  Comparison  of  figures  2  and  3  reveals  that  the 

curves  dcrivedlrom  Vai<Ay(oA)2  provide  a  more  stringent  criterion  than  those  derived  fiom  Vai(r)/(d’>^  over  the 
region  of  interest.  The  more  stringent  criterion  should  be  used  for  the  determination  of  any  tecOTmiended  operating 
conditi(Hi. 

An  additional  constraint  can  be  imposed  upon  the  variance  criteria  by  noting  that  accurate  results  are  best  obtained 
when  the  maximum  delay  time  is  approximately  four  times  the  exponential  decay  time  of  the  correlation  function 
[13].  For  the  heterodyne  case,  this  condition  reduces  to  MT At  =  4  (or  Mi^  =  4).  The  application  of  this  constraint 
permits  the  determination  of  a  single  curve  bounding  the  region  where  both  variance  ratios  are  less  than  or  equal  to 
unity.  This  curve,  which  is  presented  in  figure  4,  presents  the  value  of  r,  which  is  denoted  as  rniax(M).  that  yields 

Vai<A)/(oA)2  equaj  jq  unity  for  each  value  of  M.  For  any  specified  number  of  correlator  channels,  any  value  of  r 
less  than  or  equal  to  rinax(M)  conforms  to  the  above-mentioned  variance  criteria. 

The  curve  in  figure  4  does  not  extend  to  values  of  M  less  than  52  because  when  M  is  less  than  52  no  value  of  r 
exists  that  makes  both  variance  ratios  less  than  or  equal  to  unity.  Thus,  when  the  MT At  =  4  constraint  is  af^lied 
the  sensitivity  criterion  {ux^xised  in  this  study  cannot  be  satisfied  when  less  than  52  linear  correlator  channels  are 
used.  By  relaxing  the  acceptance  criterion  a  smaller  number  of  correlator  chaimels  can  be  accommodated.  For 
example,  if  the  maximum  acceptable  value  of  the  two  variance  ratios  is  increased  to  2  the  minimum  number  of 
linear  correlator  channels  reduces  to  M  =  26.  For  consistency,  the  maximum  acceptable  variance  ratios  are 
nevertheless  set  equal  to  unity  throughout  this  paper. 

For  applying  these  results  to  the  specific  problem  <rf  determining  acceptable  experimental  parameters  for  a  heterodyne 
DLS  experiment,  the  scattering  angle  corresptmding  to  r  =  rmaxCM), 


0max  —  2  sin  ^ 


( 


AV(fmax(M)\ 

8ii2nD  j  ’ 


(23) 


which  occurs  when  Var(A)/(oA)2  =  i,  is  plotted  for  monodisperse,  spherical,  optically  isotropic  particles  ranging  in 
size  from  30  to  90  nm  and  undergoing  a  bulk  velocity  from  0.1  to  3.5  cm/s  at  a  temperature  of  775  K  in  figure  5 
and  1300  K  in  figure  6.  This  angle  Bmax  sets  an  upper  limit  for  suitable  scattering  angles.  Equation  (23) 
incorporates  the  MTAt  =  4  constraint  mentioned  eariier.  The  relevant  flame  parameters  and  primary  particle  sizes 
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used  to  determine  the  particles'  diffusion  properties  and  the  graph's  axis  limits  were  taken  from  the  experimental 
results  of  Zhang  [15].  The  bulk  velocity  range  (0.1  to  3.5  cm/s)  was  chosen  to  be  on  the  mder  of  thennophoretic 
velocities  calculated  for  the  largest  measured  temperature  gradient  (*  228  K/mm)  of  a  6.35  mm  diameter  CO/air 
laminar  diffusion  flame  at  a  height  of  35  mm  above  the  burner  surface.  Similarly,  the  range  <rf  particle  diameters  (30 
to  90  nm)  was  chosen  to  encompass  the  size  range  of  primary  particles  observed  in  transmission  electron 
micrographs  of  iron  oxide  particles  that  were  therroophoretically  sampled  from  this  flame  when  seeded  with  Fe(CO)5 
vapor.  The  translational  diffusion  coefficient  used  to  evaluate  the  information  in  figures  5  and  6  was  evaluated  from 
the  slip-corrected  Stokes-Enstein  relation  [1, 2]; 


kT 

D  =  ^  Cs 


(24) 


where  k  is  Boltzmann's  constant,  JJ,  is  the  viscosity  of  the  solvent,  a  is  the  particle  diameter,  and 


Cs  =  l  +  -  1.234 +  0.414  exp 


(.  0,876 i)] 


(25) 


is  the  Cuimingham  correction  factor,  which  corrects  the  original  Stokes-Enstein  equatiem  for  non-continuum  effects 
when  the  gas  mean  free  path  A  is  not  negligibly  small  in  comparison  to  the  particle  diameter  a  [14].  These 
equations  indicate  that  the  slip-corrected  translational  diffusion  coefficient  increases  with  increasing  temperature. 
This  causes  0max  to  decrease  with  increasing  temperature.  This  decrease  in  Bmax  with  increasing  temperature  is 

observable  in  comparing  figures  5  and  6. 

These  Smax  results  are  qualitatively  consistent  with  resolution  considerations  for  the  heterodyne  power  spectrum 
presented  by  Ware  and  Hygare  [5],  who  define  a  frequency-based  resolution  as  the  ratio  of  the  velodty-induced 
Doppler  frequency  shift  of  the  scattered  light  to  the  diffusiem-govemed  half  width  of  the  scattered  light  spectrum; 


ii^teshift 
Aw  1/2 width  ■ 


(26) 


For  the  Lorentzian  scattered  light  frequency  spectrum  from  monodisperse  spherical  particles  [5] 


Since  this  resoluti<»  must  be  larger  than  unity  in  order  to  be  able  to  accurately  infer  the  velocity  and  diffusion 
coefficient  from  the  frequency  spectrum,  the  resolution-based  maximum  scattering  angle  has  the  following  form; 


0tnax  ~  2  sin 


(28) 


which  is  very  similar  to  the  previously  derived  expression  for  Qnnax  given  by  equation  (23). 
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4.  Application  to  Nonspherical  Particles 


The  analysis  presented  in  this  p^r  considers  the  feasibility  of  inferring  the  bulk  translational  velocities  of 
monodisperse  spherical  particles  from  heteiod)me  DLS  data.  Since  the  agglomerated  particulates  in  flame  systems 
are  often  nonspherical,  some  considerations  of  the  effects  of  ncmsphericity  and  amsotropy  to  the  uncertainties  of  the 
inferred  bulk  velocities  and  Brownian  diffusion  coefficients  are  presented. 

The  effects  of  optical  anisotropy  (and/or  nonsphericity)  on  Ii(t)  can  be  approximated  using  an  anisotropic  dipole 
model,  in  which,  for  a  cylindrically  symmetric  scatterer,  the  scatterer’s  polarizability  contains  two  independent 
components,  apar  and  ttper,  which  are  respectively  parallel  and  perpendicular  to  the  axis  of  symmetry.  The  first 
order  autocoirelation  function  Ii(t)  has  been  evaluated  for  such  scatterers  by  Beme  and  Pecora  [2]  under  the 
assumptions  that  a)  the  incident  electric  field  intensity  is  uniform  throughout  the  scattering  volume,  b)  the  position 
and  motion  of  each  scatterer  is  statistically  independent  of  the  position  and  motion  of  the  other  scatterers,  c)  the 
translational  motion  of  each  scatterer  is  statistically  independent  of  its  rotational  motion  and  vice  versa,  d)  all 
scatterers  within  the  scattering  volume  have  identical  optical  and  dynamic  properties,  e)  statistical  fluctuations  of  the 
number  of  scatterers  in  the  scattering  volume  are  negligible,  and  f)  the  length  of  each  scatterer  is  much  smaller  than 
the  wavelength  of  the  incident  light  and  thus  validates  the  use  of  the  dipole  model  instead  of  a  general  multipole 
expansion.  The  resulting  expressions  for  Ii(t)  for  W  and  VH  detection  are  [2]: 


=  {Np)  \<xf 


+  ^Pl^exp(-6DRt) 


expf-q^Dt)  exp(iq-VT)  exp(-ia)L't) 


(29) 


and 


IVh(1)  =  ejqpKq^D  +  6Dr)t:]  exp(iq  Vx)  cxp(-io)L't) 


(30) 


where  (Np)  is  the  average  number  of  scatterers  in  the  scattering  volume,  tx  =  ((Xpar  +  isotropic  part 

of  the  polarizability  tensor,  p  s  ttpar  -  otpCT  is  the  optical  anisotropy  of  the  scatterers,  and  Dr  is  the  scatterers' 
Brownian  rotational  diffusion  coefficient  The  expCiq  Vx)  term  has  been  acUedto  the  original  equations  given  by 
Beme  and  Pecora  to  account  for  the  effects  of  bulk  flow.  Likewise,  the  incident  light  angular  frequency  term  cxp(- 
which  was  neglected  by  Berafe  and  Pecora,  is  also  included  here  for  compatibility  with  the  ofiier  equations. 
Substitution  of  equations  (29)  and  (30)  into  equations  (4)  and  (6)  yields  (for  ojlo  =  wl) 
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<(n(0)n(t)Vv4»et = ^*(0))^ 


I- 


2c2<Np)  r 


NcIlO 


^  exp(-6DRt) 


+  <ii(0)>u(x) 


f  cos(q-Vt)  expC-q^Dx) 


(31) 


and 


<(n(0)n(T))vH4iet  =  <n(0))^ 


2c2{Np)P|2 
^  ISNcIiX) 


cosCqVt)  cxp[-  +  6DR)t]  +  (n(0))H(i;). 


(32) 


Since  the  W  coirelation  function  for  cylindrically  symmetric  particles  contains  two  exponentials,  the  variance 
results  of  equations  (18)  are  not  directiy  applicable  to  this  case.  Equation  (31)  demonstrates  that  assuming  a  single¬ 
exponential  form  for  the  W  correlation  function  from  cylindrically-symmetric,  nonspherical  particles  introduces  a 
systematic  error  in  the  data  interpretation.  The  percent  deviation  of  the  true  normalized  W  correlation  function,  as 
determined  through  normalization  of  equation  (31),  from  the  spherical-particle  (single-exponential)  normalized 
correlation  function,  as  given  by  equation  (9),  is  thus 


%  Error = 


F 

1  +  F 


X  100%, 


(33) 


where 


2 

exp 


andT  (s  Mit  =4/(q^D))  is  the  correlation  function's  maximum  measured  delay  time.  This  formula  indicates  that 
the  approximation  of  the  normalized  W  heterodjme  correlatitm  function  as  a  single  exponential  function  is  vahd 
when  (4lp2)/(45ict|2)  is  much  less  than  unity  and/or  (24DR)/(q2D)  is  much  greater  than  unity.  Typical  values  rf 
(4lp^)/(45lotl^)  for  prolate  ellipsoidal  particles  of  refractive  index  m  =  1.64  -  0.3i,  which  corresponds  to  FC2O3  at 
room  temperature  [16],  are  plotted  in  figure  7  for  aspect  ratios  ranging  from  2  to  20.  These  curves  were  computed 
using  the  electrostatics  j^proximations  for  the  polarizability  components  of  prolate  ellipsoids,  as  presented  by  van  ds 
Hulst  [17].  In  this  case,  the  maximum  value  of  (41p|2)/(45la|2),  which  is  obtained  as  the  aspect  ratio  approaches 
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infinity,  is  0.048.  This  limits  the  maximum  deviation  between  the  single  and  double  exponential  forms  of  the 
normalized  W  owrelationfimction  to  4.6%  for  particles  of  refractive  index  m  =  1.64  -  0.3i  at  wavelength  K  =  488 
nm.  Under  the  crmditions  in  which  the  approximation  of  the  W  heterodyne  correlation  function  as  a  single- 
exponential  function  is  deemedacceptably  accurate,  the  variance  analysis  presented  for  si^rical  particles  applies  to 
cylindrically  symmetric  particles  without  modification. 

Alternatively,  the  VH  heterodyne  correlation  function  is  a  single  exponential  function  that  has  the  same  general  form 
as  the  W  heterodyne  correlation  function  for  spherical  particles,  with  the  difference  being  the  constant  factors  and 
the  exponential  decay  rate.  The  same  general  analysis  for  9max  that  was  perf(wmed  for  the  spherical  particles  thus 

applies  for  this  case,  provided  that  F  is  replaced  with  Therefore,  when  the  constraint  MTAt  =  4  is 

used,  the  range  of  scattering  angles  tiiat  facilitate  the  inference  of  A  and  F  (s  c^D  +  6Dr)  with  variance  ratios  less 
than  or  equal  to  unity  is  found  from  the  expression: 


r 


2Tr 

A 


2iT(q^D  +  6Dr) 

q^V 


^  finax(M). 


(34) 


The  introduction  of  the  9-independent  quantity  6Dr  in  the  numerator  of  equaticm  (34)  causes  this  criterion  to  be 
quadratic  in  q.  As  a  result,  accuracy  is  not  necessarily  improved  by  running  the  experiment  at  arbitrarily  low 
scattering  angles.  When  iy<fniax(M)]^/(^*^DR)  1. where  Vq  is  the  component  of  V  in  the  direction  of  q,  two 

scattering  angles 


Snun  —  2  sin 


-4— 

[16ji^ 


\  VQrmax(M)  -  [(Vqrmax(M))^  -  96ir2DDR]i(2] 


and 


0tnax 


^  .  ll  ^ 

=  2  sin  H - :r” 


A  Varmax(M)  +  [(Vqrmax(M))^  -  96it2DDR]l/2] 


D 


(35) 


(36) 


bound  the  range  of  scattering  angles  for  which  equation  (34)  is  satisfied.  When  |y^niax(M)]^/(96nDDR)  —  1 
equation  (34)  is  only  satisfied  at  the  scattering  angle 


8  =  2  sin'^ 


A _ V(fmax(M) 

D 


16n^ 


(37) 
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When  fVifinait(M)1%96nDDR)  <  1  no  solution  exists  for  equation  (34)  and  the  variance  ratio  criterion  cannot  be 
satisfied.  The  boundaries  of  the  range  of  acceptable  scattering  angles  are  depicted  for  straight  chain  agglomerates 
suspended  in  a  gaseous  medium  at  a  temperature  of  775  K  in  figure  8  and  at  a  temperature  of  1300  K  in  figure  9. 
The  translational  and  rotational  diffusion  coefficients  used  to  calculate  the  scattering  angle  ranges  depicted  in  figures 
8  and  9  were  determined  through  the  use  of  Perrin's  equations  for  prolate  ellipsoidal  particles  of  semiaxes  a  >  b  =  c 
[4,  18,  19]: 


D  = 


kT 

6niRaj 


(38) 


and 


3kT  p2f(2-p2)G(p)-  1] 
16i^Req3  1  -  ff* 


(39) 


where 


b 


G(p)H-p=-  In 

and  Rjq  is  the  particle's  volume-equivalent  sphere  radius.  The  above  expression  for  Dr  is  valid  for  tumbling  about 
the  ellipsoidal  particle's  short  axis  b.  Assuming  that  the  particle's  material  composition  is  optically  isotropic, 
spinning  about  the  long  axis  a  does  not  affect  the  scattered  light  correlation  function.  To  account  for  non-continuum 
and  slip-flow  effects,  the  expression  for  D  is  multiplied  by  the  Cunningham  correction  factor  of  equaticMi  (25) 
(replacing  a  with  Rq).  In  the  absence  of  an  adequate  non-continuum  COTtection  factor  for  the  rotational  diffusion 
coefficient  Dr,  non-continuum  slip  flow  is  assumed  to  affect  the  rotational  diffusion  coefficient  Dr  in  a  manner 
similar  to  its  effect  on  the  translational  diffusion  coefficient  D.  Therefore,  as  in  reference  [1],  Perrin's  equation  fw 
the  rotational  diffusion  coefficient  Dr  is  corrected  for  non-continuum  and  slip  flow  ^fects  using  the  same 
Cunrungham  correction  factor  Cg  used  to  correct  the  translational  diffusion  coefficient  D. 
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Like  the  translational  diffusion  coefficient  D  for  spherical  particles,  both  the  translational  and  rotational  diffusion 
coefficients  increase  with  increasing  temperature.  As  a  result,  the  range  of  scattering  angles  conforming  to  the 
variance  criteria  decreases  with  increasing  temperature.  This  effect  is  observable  in  cmnpanng  figure  9  to  figure  8, 
where,  for  example,  increasing  the  temperature  from  775  K  to  1300  K  decreases  the  range  of  acceptable  scattering 
angles  by  36%  for  chainlike  agglomerates  of  aspect  ratio  8  undergoing  a  bulk  velocity  <rf  83  cm/s. 

As  a  result  of  the  sensitivity  considerations  of  this  section,  W  heterodyne  detection  appears  to  be  better  suited  dian 
VH  detection  for  measuring  the  bulk  velocities  of  cylindrical  agglomerates  under  the  high  temperature  conditions 
considered  here  when  the  bulk  velocities  are  less  than  approximately  30  cm/s.  For  W  detection  of  spherical  and/or 
cylindricalparticles,the  variance  ratios  of  equations  (18)  can  theoretically  be  reduced  to  less  tiian  unity  at  any  bulk 
velocity  by  reducing  the  scattering  angle  0  erf  the  detector  optics.  The  variance  ratios  from  the  VH  heteroefyne 
correlation  function,  however,  does  not  exhibit  this  behavior  because  of  the  presence  of  the  0-indcpendent  rotational 
diffusion  term  in  equation  (29).  For  VH  heterodyne  experiments  the  variance  ratios  are  both  less  than  unity  only 
within  the  limited  range  of  conditions  defined  by  equation  (31).  This  fact  makes  the  application  of  W  detection  to 
velocity  measurements  more  generally  applicable  than  VH  detection. 
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5.  Discussion 


The  expressions  developed  relating  the  variances  of  the  parameters  A  andP  of  a  heterodyne  correlation  function  to  the 
variance  of  the  noimalized  raw  data  is  applicable  to  any  discretely  samided  exponentially  decaying  cosine  function, 
{novidedthat  the  decaying  amplitude  conforms  to  a  single  exponential  function.  The  strict  validity  of  these  variance 
equations  relies  on  the  assumption  that  the  deviations  of  the  experimental  data  from  the  theoretical  form  of  g(t) 
occur  only  as  unconelated  random  errors.  Thus,  systematic  errors  are  not  considered  in  this  analysis.  In  applying 
these  results  to  the  interpretation  of  heterodyne  photon  correlation  data  from  monodisperse  spherical  particles, 
restriction  of  the  variance  ratios  of  equations  (18)  to  unity  or  less  insures  that  the  percent  uncertainties  of  the  infeired 

frequency  A  (s  q-V)  and  decay  F  (s  q^D)  terms  are  no  greaterthan  the  percent  uncertainty  (extrapolated  to  t  =  0)  cf 
the  normalized  raw  data.  This  analysis  demonstrates  that  for  accurate  data  inversion  the  ratio  r  of  the  correlation 
function's  period  to  its  decay  time  should  be  less  than  a  constant  whose  magnitude  is  typically  on  the  order  of 
magnitude  of  10.  At  larger  values  of  this  ratio  the  correlation  function  decays  to  its  baseline  value  before  it  can 
complete  enough  of  an  oscillation  to  contain  enough  information  for  accurate  frequency  inference. 

These  variance  results  can  be  used  to  estimate  the  uncertainties  of  Vq  and  D  through  the  relationships  (see  [20]): 


and 


V»(D).(^]  ''v.KD*(f] 


Vai<q). 


(40) 


(41) 


The  statistical  uncertainties  of  VqandD  are,  by  definition,  the  square  roots  of  their  variances.  The  variances  of  A 
and  r  requited  for  evaluating  equations  (40)  and  (41)  are  readily  obtained  ftom  equations  (18).  The  variance  of  q, 
which  increases  with  decreasing  scattering  angle,  is  determined  by  both  the  angular  resolution  of  the  goniometer  used 
to  set  the  scattering  angle  and  the  size  of  the  scattering  volume  [7]. 

As  previously  discussed,  the  variance  expressions  presented  in  this  paper  can  also  be  applied  to  cylindrically 
symmetric  particles.  Since,  unlike  the  VH  heterodyne  correlation  function,  the  W  correlation  function  ftom 
monodisperse  cylindrical  particles  contains  a  multi-exponential  decay  factor,  equations  (18)  are  not  strictly  applicable 
since  the  approximation  of  the  correlation  functicHi  as  a  single-exponential  function  introduces  a  systematic  error  in 
the  data  analysis.  The  correlation  function  does,  however,  diten  follow  the  single-exponential  form  very  closely 
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under  the  conditions  presented  in  the  previous  section.  One  possiWe  extension  of  the  present  analysis  would  be  to 
ifetermine  how  this  systematic  error  affects  the  variance  estimates  of  the  inferred  parameters. 

The  riinctional  form  of  the  correlation  function  also  deviates  from  the  single-exponential  form  of  equation  (9)  when 
the  scattering  particles  are  polydisperse,  thus  exhibiting  variations  in  the  particle  dimensions  and/or  optical 
properties.  Like  the  case  of  W  heterodyne  detection  from  monodisperse,  nonspherical  scatterers,  this  deviation 
produces  a  systematic  error  in  the  data  analysis  if  a  single-exponential  decay  is  assumed.  The  theoretical  correlation 
function  of  light  scattered  from  a  suspension  of  polydisperse  scatterers  is  generally  expressed  as  a  weighted  integral 
of  the  monodisperse  correlation  function  over  the  entire  range  of  particle  dimensions  and/or  (q)tical  properties.  In 
order  to  properly  evaluate  the  effects  of  polydispersity  on  the  statistical  uncertainty  of  the  inferred  results,  an 
inversion  scheme  that  accounts  for  this  polydispersity  must  first  be  developed  for  the  heterodyne  case.  Then,  a 
variance  analysis  could  be  performed  using  the  same  general  methodology  as  the  analysis  presented  in  this  paper. 

One  further  limitation  of  this  analysis  is  the  assumption  that  the  pre-exponential  factor  of  equation  (8),  (31),  or  (32) 
is  known  prior  to  the  inversion  of  the  raw  data  so  that  the  data  can  be  properly  normalized.  This  factor  can  be 
determined  by  extr^x)lating  the  correlation  function  to  its  2erx><lelay-time  value.  If  the  current  formulation  is  used  to 
estimate  the  variances  of  A  andF  from  the  variances  of  the  correlation  function  data,  the  variance  of  the  normalized 
correlation  function  must  first  be  adjusted  to  account  for  the  uncertainty  of  this  pre-exponential  factor.  An 
alternative  technique  would  be  to  treat  the  pre-exponential  factor  as  a  third  unknown  parameter  in  the  least-squares 
data  inversion  scheme  and  extending  the  current  variance  analysis  using  a  derivation  similar  to  the  one  presented  in 
this  paper. 
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6.  Summary  of  Results 


The  results  of  the  present  study  may  be  summarized  as  follows: 

a)  As  the  particle  size  and  bulk  velocity  increases,  the  maximum  ^x^eptable  scattering  angle  for  accurate  determination 
of  the  bulk  velocity  and  Brownian  diffusion  coefficient  increases.  This  occurs  because  increasing  the  size  or  bulk 
velocity  increases  the  decay  time  of  the  correlation  function  with  respect  to  its  period,  thus  permitting  more 
measurable  oscillations  of  the  cosine  factor  before  the  correlation  function  completely  decays  to  its  baseline  value. 

b)  On  the  basis  of  the  variance  analysis  for  monodisperse  sjAerical  particles,  the  inference  of  A  (s  q-V)  and  F  (=  q^D) 
with  a  percent  uncertainty  less  than  or  equal  to  the  initial  percent  uncertainty  of  the  normalized  COTelation  function 
data  may  be  performed  for  the  CO/air  diffusion  flame  under  consideration  when  the  average  particle  velocity  is  greater 
than  approximately  1  cm/s  and  the  diameter  of  the  particles  is  larger  than  approximately  40  nm.  The  translational 
diffusion  coefficient  D  and  the  q-component  of  the  average  velocity  Vq  may  then  be  found  through  knowledge  of  q 
from  equation  (2). 

c)  Under  certain  conditions,  this  analysis  can  be  extended  to  non~spherical  particles.  For  W  detection  of  monodisperse 
cylindrically  symmetric  particles,  the  heterodyne  correlation  function  contains  two  distinct  exponential  terms, 
making  the  single-exponential-based  variance  analysis  consideredin  this  paper  invalid.  When  the  scattering  particles’ 
polarizability  anisotropy  is  much  smaller  than  the  isotropic  part  of  the  polarizability  tensor  and/or  the  rotational 
relaxation  time  is  much  smaller  than  the  translational  relaxation  time,  the  W  correlation  function  may  be 
approximated  as  a  single-exponential  function.  In  this  case,  the  variance  analysis  presented  in  this  paper  applies  to 
non-spherical  particles  without  modification.  The  variance  analysis  also  applies  to  the  single-exponential  VH 
correlation  function  from  cylindrically  symmetric,  anisotropic  particles.  In  this  case,  however,  the  addition  of  a 
scattering-angle-independent  rotational  diffusion  term  to  the  exponential  decay  factor  changes  the  qualitative 
characteristics  of  the  range  of  scattering  angles  giving  acceptable  error  propagation  characteristics.  For  particles  of 
aspect  ratio  less  than  approximately  ten  the  magnitudes  of  the  minimum  measurable  bulk  velocities  are  ccMisiderably 
larger  (by  a  factor  of  10  to  100)  than  the  velocities  measurable  using  W  detection.  As  revealed  in  figure  8  the 
minimum  measurable  bulk  velocity  gradually  deoeases  as  the  aspect  ratio  of  the  scatterers  increases.  This 
occurrence  of  a  distinct  minimum  measurable  bulk  velocity  severely  limits  the  ^plicability  of  VH  heterodyne 
detection  to  measuring  bulk  velocities. 
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Figure  1 ;  Scattering  geometry  for  the  dynamic  light  scattering  (DLS)  expserimenL 


No.  of 
Correlator 
Channels 


Figure  2:  Boundaries  of  the  acceptable  operation  region,  corresponding  to  Var(A)/(aA)2  -  1,  for 
the  inference  of  A  (■  q-V)  from  a  single-exponential  heterodyne  DLS  correlation 
function  determined  for  M  linear  correlator  channels.  The  system  should  operate  below 
these  curves. 
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Period/(Decay  Time),  (25rr)/A 


No.  of 
Correlator 
Channels 


(Sample  Time)/Period,  (AAt)/(27r) 

Figure  3;  Boundaries  of  the  acceptable  operation  region,  corresponding  to  Var(r)/(ar)2  =  1,  for 

the  inference  of  F  (■  q^D)  from  a  single-exponential  heterodyne  DLS  correlation 
function  determined  for  M  linear  correlator  channels.  The  system  should  operate  below 
these  curves. 
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Figure  7-  Relative  optical  anisotropy  parameter  (4ipp)/(45|ap)  of  prolate  ellipsoidal  particles  of 
refractive  index  m  *  1.64  -  0.3 i  (corresponding  to  room  temperature  FeiOs)  as  a 
function  of  aspect  ratio. 
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VII.  Computer  Source  Code  (  AGGL  )  Developed  for  the  Prediction  of  the 
Light  Scattering  Properties  of  Agglomerated  Particulates  Composed  of 
Rayleigh  Sized  Primary  Particles. 

The  FORTRAN  codes  inciiwifd  here  evolved  from  the  Exact  Formulation  for  the  light  scattering  properties  of 
agglomerated  particulates.  The  exact  formulation  as  noted  in  section  I  of  this  Report  has  been  published  in  the 
Journal  of  Physics  D:  AppUed  Physics  Vol.  27,  pp  2258  -  2270(1994).  A  Guide  about  (he  running  of  this  code  and 
an  output  example  are  provided  in  section  VI  of  this  Report  More  information  about  this  code  may  be  found  in  the 
Masters  Thesis  of  Wujiang  Lou  entiUed  ‘Theory  and  Application  of  eectromagnctic  Wave  Scattering  and 
Absorption  of  Agglomerated  Small  Spherical  Particles”  .  This  document  is  available  in  the  Louisiana  State 
University  Library. 


c.a.  sub:  a  collection  of  all  subroutines  ‘em 

c.a.  ‘e™ 

C 

C23456789012345678901234567890123456789012345678901234567890123456789 


c^*****  preparedby  w.  lou,  dccember  1994 
Q —  version:  d- featuring  derivatives  of  scattering  quantities 


c — sub:  aggLshell 

Q —  a  pilot  program  providing  interface  to  master  program 

c - — - lou,august  10,1994 

subroutine  aggl_shell(np,ns,nps,alpha,rfindx) 

parameter  (  mnp=1004id=181) 

complex  rf!ndx,eps]on 

diaiacter  eff*  1  ,poI*  1  ,chain*  1 

real  kvve(nd),  khlie(nd)Jcvbe(nd)4dive(nd), 

$  qetn,qetk,dkvvn(nd),dkvvk(nd),dkhvn(ndXdkhvk(nd), 

$  qehtn,qehtk,dkhlm(nd),dldihk(nd),dkvhn^^ 

$  qetx,qehtx,dkvvx(nd),dkhvx(i]d),dkvhx(nd),dkhhx(nd) 
common/p2/pol,eff, chain 

common/j2s/scahh(181),scavv(181),  scahv(181X  scavh(181), 
$  cosang(181),sinang(181) 
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common/l  l/xo(mnp),yo(nmp),zo(mnp) 
conunon/Iou/lu^u  1  »lu2,Iu3  4u4 

comraon/loul/cospsi(10)^inpsi(10),coschi(10),sinchi(10),w(10) 

common/orient/dchi,(4>si 

coinnion/q(jext/qcxtO,qabsO,qscaO,qexthO,qabdiO,qscahO 
commoii/loa/qct2iO,qctlrf),dkvvn(Xnd),dkvvW)(nd),dto 
$  dkhvk(XndXqehtiiO,qchtkO,dkhhn{Knd),dk^^ 

$  dkvhiiO(Dd),dkvhkO(nd) 
cominon/lox/qetxO,dkvvxO(nd),dkhvxO(nd), 

&  qchtxO,dkhhxO(ndXdkvlu^ 

Q - averaged  derivatives  for  random  structures/pass  to  master 

common/qextm/qext,qabs,qsca,qexth,qabslM}scah 
common/j2sm/scahhk(181),scavvk(180^  scahvk(181),  scavhk(181), 
$  angk(181) 

common/iouv/rvv(Dd),drvvn(nd),drs^k(nd),dkhvcn(ndXdkhvdc(nd), 

$  kvve,dkvven(nd),dkvvek(nd)^v(ndXditivn(ndX 
$  driivk(nd)Jdive,kvhe 

common/Iouh/ihh(nd),drfilm(nd),drf^  dkvhcn(ndXdkvhek(nd), 
$  khhe,dkhhen(ndXdkhhek(nd)4vh(iidXdrvhn(ndXdr^ 

common/Ioux/qctx,drvvx(nd),drlivx(ndX  drvrhx(nd)4rfibx(ndX 

$  qehtx,dkvvex(nd),dkhvcx(nd),dk\iicx(nd),dldiliex(nd) 

if(chain.eq.V)then 

else 

nps=l 

end  if 

pi=3.1415926 

cpslon=ifindx*ifindx 

c — data  initialization - - 

call  loul 

sumxl=0.0 

qsca=O.0 

qcxt=0.0 

qabs=O.0 
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qscali=0.0 
qextb=0.0 
qabsh=0.0 
do  96  i=l,  ns 
scawk(i)=0.0 
scahhk(i)=0.0 
scavlik(i)=0.0 
8cahvk0)=0.0 
anglc=angk(i)*pi/180.0 
cosangO)=cos(ang]e) 
sinang(i)=sin(angle) 

96  continue 
qetn=0.0 
qetk=0.0 
qehtn=0.0 
qehtk=0.0 
qete=0.0 
qehtx^.O 
do  97  i=l,ns 

dkwn(i)=0.0 

dkwk(i)=0.0 

dkhvn(i)=0.0 

dkhvk(i)=0.0 

dktihn(i)=0.0 

dklihk(i)=0.0 

dkvhn(i)=0.0 

dkvlikCi)=0.0 

dkvvx(i)=0.0 

dldivx(i)=0.0 

dkvhx(i)=0.0 

dki]hx(i)=0.0 

97  continue 

do  100  kk=l,  nps 
if  (chain  .eq.  V)  then 
call  ran<kHn(np) 


else  if  (chain  .eq.  's')  then 
call  strait(np)- 
else  if  (chain  .eq.  'c')  then 
call  clustKnp) 
dse 

print  *,  'no  such  choice  of  chain' 
stop 
endif 

if(lu3.eq.0)then 

if(iu4.eq.0)then 

call  agglou(Knp^pha,epslon^) 
else 

call  agglou0d(np,alpha,ifindx4is) 
endif 
else 

chi=dchi*pi/180.0 

psi=dpsi*pi/180.0 

call  aggloul(np^pha^ndx4is,chi,psi) 
endif 

qsca=qsca(Vnps-t-  qsca 
qext=qext+  qextO/nps 
qabs=qabs+  qabsO/nps 
qscah=qscahO/nps+  qscah 
qexth=qexth+  qexthO/nps 
qabsh=qabsh+  qabdiO/nps 
sumxl  =  sumxl  +  agxl/nps 
do  98  nini=l^ 

scavvk(mm)=scawk(inm)+scavv(inni)/nps 

scahhk(mm)=scahhk(inm)+scahh(iran)/nps 

scavhk(nim)=scavhk(inm}+scavh(mm)/nps 

scahvk(nnn)=scahvk(niin)+scahv(nini)/nps 

98  continue 

qetD=qetn-fqetn(Vnps 

qetk=qetk4qetk0/nps 


qehta=qehtiH^htaQ/nps 

qelitk=qehtk-H]ehtkO/nps 

q^=qetx-fqetx(Vnps 

qehtx=qehtx+qebtxO/nps 

do99i=l4»s 

dlcvvnO)==dkvvn0)+dkvvnOO)/nps 

dkvvk(i)==dkwk0)+dkvvkOO)/nps 

dkhvnG)=dkhvn(i)+dkhvnO0)/nps 

dkhvk(i)=dkh^*O>+dkhvkO0)/nps 

dkh}mO)=dldtbn0>+<lldihiiO(i)/nps 

dkhhk(i)=didM0>f<lkiM()(iyiips 

dkvhii(i)=dkvhnQ)+dkvhiiO(i)/nps 

dkvhk(i)=dkvhk(i)+dkvhkO(i)/nps 

dkvvx0)=dkvvx0)+dkvvx()0)/nps 

dkhvx0)=dkhvxQ>KlkhvxO(i)/nps 

dkvhxO)=dkvhxO)+<llcvhxO(i)/nps 

dkhhx(i)==dkhlixQ)+dkhhxO(i)/nps 

99  continue 

100  continue 

c — derivatives - 

do  211  k=l^s 
kvve(k)=scavvk(k)/qext 
kbve(k>=scabvk(kyqext 
dfcvven(k)=(dkvvn(k)-kvve(k)*qetnyqext 
dkvvek(k)=(dkv^4(k>kvve(k)*qeflcyqext 
dkhven(k)=(dkhvn(k>-khve(k)*qeuiyqext 
dkhvek(k)=(<^vk(k)-khve(k)*qetkyqext 
dfcvvex(k)=(dkvvx(k)-kvve(k)*qeftyqcxt 
dkhvex(k)=(dkhvx(k)-khve(k)*qelxyqext 
rw(k)=scavvk(kyscavvk(ns-k+l) 
iliv(k)=scahvk(kyscavvk(k) 

211  continue 

do  230  k=:l^ 

drvvn(k)=(dkvvn(k)-rvv(k)*dkvvn(ns-k+l)yscavvk(ns-k+l) 
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d^vvk(kK<^kv\^)^Tv(k)*dkv^1c(m-k+l)yscaw 

dii»vn(kMdkhvn(k)-iiiv(k)*dkvvii(k)yscavvk(k) 

dilivk(k)=(dkhvk(k)-ihv(k)*dkvvk(k)yscavvk(k) 

<lrvvx(k)=(<fl'^^(k)-rvv(k)*dkvvx(ns-k+l)yscav>ic(ns-k+l) 

driivx(k)=(<lkhvxO'^)-Av(k)*<Ikvvx(k)yscav\ic(k) 

230  continue 
do  311  k=Ms 
klihe(k>=scahhk(kyqextfa 
kviie(k)=scavhk(kyqextii 

dkhhek(kKdkhW^)-l*^ 

dkvhen(k)=(<ikvl“>(k>4vhe(k)*qehlnyqexth 

dkvhek(k)=(dkvhk(k>4cvhe(k)*qehtkyqexfh 

dkhhex(kM<ffl^W«O^Hi^e(k)*qehtxyqexth 

dkvIiex(k)=<<JkvlK(k)-fcvhe(k)*qehlxyqexth 

riih(k)=8calihk(kyscahhk(ns-k+l) 

rvh(k)=scavhk(kyscahhk(k) 

311  continue 

do  330  k=:l4)s 

diMui(kMtikhhn(k)-ilifa(k)*dkhlm(ns-k+l)yscahhk(ns-k+l) 

driihk(kMdklM{kyitJi(k)*dkhhk(ns-k+l)yscahhk(ns-k+l) 

drvlm(k)=(<fl^vbnO^)'*vh(k)*dkhhn(k)yscahhkOt) 

drvlik(k)=(<lkvhk(k)-rvh(k)*dkhhk(k)yscahlik(k) 

drijhx(k)=(dkhhx(k>-riih(k)*dkhlix(ns-k+l)yscahhk(ns-k+l) 

divhx(k)=<<Uwhx(k)-rvh(k)*dkIihx(k)yscahhk(k) 

330  continue 
letiun 
end 


c- - - - — - 

c — sub:  aggloul 

c - ^forchains  at  a  specific  orientation 

c - dateof  generation: July  15,  1994  —by  wujiang lou 

c - derivatives  included  in  august  1994  - 
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c:  qetn(kXx):  derivatives  of  qext  wrt  n,k,  and  xp  respectively 
c:  qehtn(kXx):  same  as  qetn(kXx)  but  for  vertical  incidence 
c:  dkabn(kXx):  derivatives  of  differential  scattering  efficiency 
c -  wrt  njc^d  xp.  a&b  denote  polarization(v  or  h) 

subroutine  agglou  l(np,a]pha/findx4is,chi,psi) 
panuneter  (mnp=l(X)4id=181) 

complex  rfindx,epslon,e(34iuip),esum/sum,ep(6),theta, phi, 

$  cm,  eh(3,mnp),  eph(6),  esumh4'sumh,thetah,phih,  ci, 

$  el(3,mnp),ehl(3,mnp),ex(3,mnp),ehx(3,mnp),  ctem 
character  pol*l,chain*  1  ,eff*  1 
common/p2/pol,eff,chain 
common/lou/lu,lul,iu2,lu3,lu4 
common/jl/e,  di 
common/ji2/el,  ehl 
common/jI3/ex,ehx 
common/12/x(mnp),y(mnp),z(mnp) 
common/ll/xO(mnp),yO(mnp),zO(mnp) 
common/j2s/scalih(181),  scavv(181),  scahv(181),  scavh(181), 

$  cosang(181),sinang(181) 
common/qqext/qext,qabs,qsca,qexth,qabsh,qscah 
cornmon/loa/qetn,qetk,dkvvn(nd),dkvvk(nd),dkhvn(nd),dkhvk(nd), 
$  qehtn,qehtk,dkhhn(nd),dkhhk(nd),dkvhn(nd)4kvhk(nd) 
common/lox/qctx,dkvvx(nd),dkhvx(nd),qehtx,dkhhx(nd),dkvlK(nd) 

data  py3. 1415926/ 

epslon=rt!ndx*rtindx 

ii]{dia2=a]pha*alpha 

ci=cmplx(0.0,1.0) 

if(lu.eq.2)then 

yjl=(sin(alpha)/alpha-cos(alpha))/alpha 

else 

yjl=alpha/3.0 
end  if 

call  cootd(np4dpha,l.,0.0,l.,0.0) 
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call  zhoul(np,aIpha,cpsIon) 


chicos=cos(chi) 

chisin=sin(chi) 

psicos=cos(psi) 

psisiD=sin(psi) 

call  coord(np,alpIia,clucos,chisiii,psicos,psism) 
if(cfT.cq.*n*  .or.  cff.cq.'iOtben 
call  zhou0(np,alpha,epslon) 
else 

if(lii4.eq.0)then 

call  zhou3(np,alpha,epslon) 
else 

call  zhou31(np^pha4findx) 
end  if 
end  if 

c — extinction  &  absorption  effiency  factors 
asuin=0.0 

esiiin=(0.,0.) 

asuinh=0.0 

csumh=(0.,0.) 

do40j=l,np 

ctem=cmplx(cos(z(j))rsm(z(j))) 
csum=esum+  c(l  j)*ctem 

asiun=asuiD+e(lj)*conjg(e(lj))+«(2j)*conjg(e(2j)) 

$  4c(3j)*conjg(e(3J)) 

esuiiih=:csumh+eh(2j)*cteni 

asunfih=asunih+eh(  1  j)*conjg(eh(  1  j))+ch(2  j)*conjg(eh(2  j)) 
$  +eh(3j)*conjg(eh(3j)) 

40  continue 

qext;=aiinag((epsIon-l  .)*esum) 
qabs=aimag(epslon- 1 .0)*asuin 
qcxtb=aimag((epsIon“l.)*csuinh) 
qabsh=^mag(epsIon- 1 .0)*asumh 
qext=4.*yjl/np*qcxt 
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qabs=4.*yjl/np*qabs 
qextb=  4.*yj  l/np*qexth 
qabsh=4.*yj  l/np*qabsh 

e — ^total  scattering  cffiency  factor,  set  lul=0  to  skip 
if(1ul.eq.0)thea 
qsc3=qext-qabs 
qscali=qext])-qabsh 
else 

call  louO(np,alpha,epslon,qsca,qscah) 
end  if 

c — angular  differential  scattering  eiTidency 
do  140  k=l,ns 
cangle=cosang(k) 

8angle=sinang(k) 
do  120 1=1,3 
ep(l)=cmplx(0.,0.) 

120  cph(l)=cinplx(0.,0.) 
do  132  j=l,np 
rj=z(j)*cangle+y(j)*sangle 
ctem=cniplx(cos(ii),-sin(ij)) 
do  130 1=U 
ep(l)=ep(l>+ctein*eOj) 

130  q>h(l)=eph(l)+ctein*eli(lj) 

132  continue 

theta=cangle*ep(2>sangle*ep(3) 

phi=ep{l) 

scahv(k)=theta*conjg(tfaeta) 

scaw(k)=plii*conjg(phi) 

thetali=cangle*eph(2)-sangle*eidi(3) 

phili=eph(l) 

scalih(k)=thetah*conjg(dietah) 

scavh(k)=pIiih*conjg(pliih) 


140  continue 


c=(cabs(epslon-l.))**2*alpha**2*yjl*yjl/pi/np 

do  170  k=Ms  ■ 

scahh(k)=c’''scahh(k) 

scavh(k>=c*scavh(k) 

scahvOc)=c*scahv(k) 

170  scaw(k)=c*scaw(k) 

c — partial  derivatives  of  extinction  and  differential  scattering 
if(lu4.eq.0)tfaen 
return 
end  if 

c - vertical  incidence:  qext,cvv  and  chv  wit  (n,k,xp=alpha) 

csuin=(0.0,0-0) 

fsum=(0.0,0.0) 

do250j=l4»P 

cin=cniplx(cos(z(j)),-sin(z(j))) 

esuin=esum+  (e(  1  J)*2.0*rfindx+(epslon- 1 .0)*el  (1  j))*cni 

fsuin=fsuin+(ex(lj)-cmpIx(O.OXj))*e(ljValpha)*cm 

250  continue 

qetn=4.0/np*yj  1  ♦aiinag(esuin) 
qetk=4.0/np*yj  1  *ieal(esum) 

qetx=(1.0-0.3*alpha2)/alpha/(1.0^.1*alpha2)*qext 
&  +  4.*alpha/3.0/np*aiinag((epslon-1.0)*fsum) 

tcinp0=(epslon-1.0)*conjg((epslon-1.0)) 
templ=alplia2*yjl*yjl/pi/np 

teinim==4.0*real(rfmdx)*(rfindx*conjg(ifiiidx)-1.0) 

tempk=4.0*ainiag(ifindx)*(rfindx*conjg(rfindx)+1.0) 

do  210  k=l4is 

cangle=cosang(k) 

sang]e=sinang(k) 

esuin=(0.0,0.0) 

fsuni=(0.0,0.0) 

ep(l)=(0.0,0.0) 

ep(2)=(0.0,0.0) 

ep(5)=(0.0,0.0) 
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eph(5M0.0,0.0) 
do220j=l4ip  ■ 

ij=z(j)*cangIe+yO)*sangle 

cin=cmpIx(cos(ij),-sin(q)) 

esuni=esum-(-  c(l  j)*cin 

fsuiii=fsum+  el(l  j)*cin 

cp(l)=ep(l>t<cangle*c(2j)-sangle*e(3j))*cm 

ep(2>=ep(2>+(cangle*el(2j)-sangle*el(3J))*cm 

ep(5)=ep(5>+<ex(lj)<i*ij/alpha*e(lj))*cm 

eph(5)=eph(5>+(cangle*(ex(2j>e(2j)*ci*ij/alpha) 

&  -sangle*(ex(3j)-c(3j)*ci*ij/alpha))*cm 

220  continue 

thefa=esum*conjg(fsum)+conjg(esum)*fsuni 
cp(3)=ep(l)*conjg(epC2)>+conjg(ep(l))*ep(2) 
phi=cniplx(0.0,10)*(conjg(esuin)*fsuna-esuni*conjg(fsuin)) 
ep(4)=cmplx(0.0,1.0)*(conjg(ep(l))*ep(2)-cp(l)*conjg(ep(2))) 
dkvvn(k)=scavv(k)/teinpO*tempn+tempO*templ  ♦theta 
dkwk(k)=scaw(k)/tempO*tempk+tempO*templ*phi 
dkhvn(k>=scahv(k)/tenipO*tempn+tempO*templ  ♦ep(3) 
dkhvk(k)=scahv(k)/tempO*tenipk+tempO*tcmpl  *ep(4) 
dkvvx(k)=(4.0-1.2*alpha2)/alplia/(l-0.2*alpha2)*scaw(k) 

&  +tempO*templ  *(ep(5)*conjg(esum)+conjg(ep(5))*esum) 

dkhvx(k)=(4.0-1.2*alpha2)/alpha/(l-0^*alpha2)*scahv(k) 

&  +tempO*templ  *(eph(5)*conjg(ep(l)>+conjg(eph(5))*ep(l)) 
210  continue 

c— derivatives  of  qexth,  cvh,  &  chh:  horizontal  incidence. 
esunih=(0.0,0.0) 
fsumh=(0.0,0.0) 
do350j=Mp 

cni=cniplx(cos(z(j)),-sin(z(j))) 

esunih=esunib+  (eh(2  j)*2.0*ifindx4<cpslon- 1 .0)*eh  1  (2  J))*cm 

fsuinh=fsunih+(ehx(lj)-cmplx(O.OXj))*eh(ljyalpha)*cm 

350  continue 

qehtn=4.0/np*yj  1  *ainiag(esumh) 
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qelilk=4.0/np*yj  1  *rea!(esumh) 

qeta=(1.0-0.3*alpha2)/alpha/(1.0-0.1*alpha2)»qexth 
&  +  4.*alplia/3.0/np*aiinag((cpslon-1.0)*fsunih) 

: - ^angular  derivati  ves- 

do330k=l,ns 
caiigIe=cosang(k) 
sangle=sinang(k) 
esuinh=(0.0,0.0) 
fsuiiih=(0-0»0.0) 
eph(l)=0.0 
ci*(2)=0.0 
ep(5)=(0.0,0.0) 
ei*(5)=(0-0A0) 
do320j=Mp 
q=z(j)*cangle+y(j)*sangle 
cin=cmplx(cos(g),-sin(ij)) 
esuinh=esuinh+  eh(l  j)*cni 
fsun[ib=fsuinh+  ehl(l  J)*cm 
eph(l)=eph(l>Kcangle*eh(2j)-sangle*eh(3j))*cni 
eph(2)=eph(2}+<cangle*ehl(2j)-sangle*chl(3j))*cni 
ep(5)=ep(5>Keltt(lj)-ci*i37alpha*ch(lj))*cm 
q)h(5)=eph(5)+<can8Je*(ehx(2j)-eh(2j)*ci*ij/aIpha) 

&  -sangle*(ehx(3j)-eh(3j)*ci*ij/alpha))*cni 

320  continue 

lhetah=esuinh*conjg(fsunjh)+conjg(esunih)*fsuinh 
ei*(3)=eph(l)*conjg(eph(2)>+conjg(eph(l))*eph(2) 
phUi=cmplx(0.0,1.0)*(c(Mijg(esuinh)*fsuinh-esuinh*conjg(fsuinh)) 
eph(4)=cmplx(0.0,l  .0) 

$  *(conj^eph(l))*eph(2)-eph(I)*conjg(eph(2))) 

dkvhnGc)=scavh(k)/tempO*tempn+tempO*templ*thetah 
dk>1ik(k)=scavh(k)/temp0*tempk+temp0*templ  ♦phih 
dkhhn(k)=scahh(k)/tenip0*tempn4<enip0*templ  *cph(3) 
dkhlik(k)=scahh(k)/tenip0*tenq)k4tenip0*tenipl*cph(4) 
dkvlix(k)=(4.0-1.2*alpha2yalpha/(l-0.2*alpha2)*scavh(k) 

&  +terapO*templ  *(cp(5)*conjg(esunih>Konjg(ep(5))*esumh) 
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dkhhx(k)=(4.0-1.2*alpha2)/alpha/(l-0.2*alpha2)*scahh(k) 

&  +temp0*templ*(eph(5)*conjg(eph(l)>+conjg(eph(5))*eph(l)) 
330  continue 
return 
end 
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c¬ 


c — sub:agglouOd 

c - agglouO  with  derivatives 

c - dateof  generation:  august  20,  1994 

c - - — - - — - -  - 

subroutine  agglou0d(np,alpha,rfindx4is) 
parameter  (mnp=100,nd=181) 

complex  rfindx,epslon,e(34nnp),esum/sum,ep(6),theta,phi, 

$  cm,  eh(3,nmp),  eph(6),  esumh,fsumh,thetah,phih,ci, 

$  el(3,mnp),chl(3,mnp),ex(3,mnp),ehx(34nnp),  ctem 
character  pol*  1  ,chain*  1  ,eff*  1 
real  scahh(nd),scavv(nd),scahv(nd),scavh(nd), 

$  qeta,qetk,dkvvn(nd),dkvvk(nd),dkhvn(nd),dkhvk(nd), 

$  qehtn,qehtk,dkhhn(nd),dkhhk(nd),dkvhn(nd),dkvhk(nd), 

$  qetx,qehtx,dkvvx(nd),dkhvx(nd),dkvhx(nd),dkhhx(nd) 
common/p2/pol,eff,chain 
common/lou/lu,luUu2,lu3,lu4 

common/loul/cospsi(20),sinpsi(20),coschi(20),sinchi(20),w(20),ng 

common/jl/e,  eh 
common/jl2/el,  ehl 
common/jl3/ex,ehx 
common/12/x(mnp),y(mnp),z(mnp) 
common/1  l/x0(mnp),y0(mnp),z0(mnp) 
c — orientation  averaged  quantities 

common/j2s/scahh0(181),  scavv0(181),  scahv0(181),  scavh0(181), 
$  cosang(181),  sinang(181) 

common/qqext/qextO,qabsO,qscaO,qexthO,qabshO,qscahO 

c(Mnmon/Ioa/qetnO,qedcO,dkvvnO(nd),dkvvkO(nd),dkhvnO(nd), 

$  dkhvkO(nd),qehtnO,qehtkO,dkhhnO(nd),dkhhlcO(nd), 

$  dkvhn0(nd)4kvhk0(nd) 
common/lox/qetxO,dkvvxO(nd),dkhvxO(nd), 

$  qehtxO,dkhhxO(nd),dkvhxO(nd)  . 

pi  =  3.1415926 
pi2=pi*0.5 
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alpha2=alpha*alplia 

epslon=ifindx*ifindx 

ci=cinplx(0.0,1.0) 

if(lu.eq.2)then 

yj  l=:(sin(alpha)/alpha-cos(aipha))/aIpha 
dse 

yjl=alpha/3.0 
end  if 

call  coord(np^pha,l.,0.0.1-i0.0) 

call  zhoul(Dp,alpha,epsIon) 

qext0=0. 

qsca0=0. 

qabs(>=0. 

qextli0=0. 

qscah0=0. 

qabsh0=0. 

qetn0=0.0 

qetk0=0.0 

qehtD0=0.0 

qehtk0=0.0 

c — data  initialization - 

do  10k=l,ns 
scalib0(k)=0. 
scahv0(k)=0. 
scaw0(k)=0. 
scavh0(k)=0. 
dkvvn0(k)=0.0 
dkwk0(k)=0.0 
dkbvn0(k)=0.0 
dkh\ic0(k)=0.0 
dkbhn0(k)=0.0 
dkhlik0(k>=0.0 
dkvlin0(k)=0.0 
dkvlik0(k)=0.0 


10  continue 


do420ip=l4ig 

I)sicos=cospsi(ip) 

psisin=sinpsi(ip) 

do410ic=Mg 

chisin=sinchi(ic) 

chicos=coschi0c) 

call  coord(np,alpha,chicos,chisin,psicos,psisin) 

if(eff.eq.'n'  .or.  eflF.eq.V)then 
call  zhou0(np,alpha,epslon) 
else 

call  zhou31(np^pha/findx) 
end  if 

c — extinction  &  absorption  eftiency  factors 
asuni=0.0 

esum=(0.,0.) 

asuinh=0.0 
esunih=(0.,0.) 
do40  j=l4»P 

cteni=cm{rfx(cos(z(j)),-sin(z(j))) 
esum=esuni+  e(l  j)*ctem 

asuin=asuni4e(l  j)*conjg(e(l  j))+€{2 j)*conjg(e(2  j)) 

$  +e(3j)*conjg(c(3j)) 

esumh=esumh+ch(2j)*ctem 

asunih=asunih+eh(lj)*conjg(eh(lj)>+ch(2j)*conjg(eh(2j)) 

$  4eh(3j)*conjg(eh(3J)) 

40  continue 

qcxt=aiinag((cpslon-l.)*esum) 

qabs=aimag(epslon-1.0)*asum 

qexth=ainiag((epsIon-l.)*esuinb) 

qabsh=ain]ag(epslon-1.0)*asumb 

qext=  4.*yj  l/np*qext 

qabs=4.*yjl/np*qabs 

qexth=  4.*yj  l/np*qexth 
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qabsh=  4.*yj  l/np*qabsh 

c — ^total  scattering  effiency  factor,  set  lul=0  to  skip 
ifGul.eq.O)then 
qsc3=qext-qabs 
qscah=qexth-qabsh 
else 

call  IouO(Dp,alpha,epslon,qsca,qscah) 
end  if 

c — angular  differential  scattering  efiidency 
100  continue 
do  140  k=l,ns 
cangle  =  cosang(k) 
sangle  =  sinang(k) 
do  120 1=1,3 
epO)=cniplx(0.,0.) 

120  eph(l)=cmplx(0.,0.) 

do  132j=l,np 
rj=z(j)*cangle+y(j)*sangle 
ctein=cmiHx(cos(i3),-sin(ij)) 
do  130 1=U 
cp(l)=ep(l>+«tein*e(lj) 

130  ephO)=cph(l>+ctem*ehOj) 

132  continue 

theta=cangle*ep(2)-sangle*ep(3) 

plii=ep(l) 

scahv(k)=tlieta*conjg(theta) 

scaw(k)=phi*conjg(phi) 

thetah=cangle*eph(2)-sangle*eph(3) 

pbih=eph(l) 

scahh(k)=tbetah*conjg(thetali) . 
scavh(k)=phih*conjg(phih) 

140  continue 

c=(cabs(epslon-l.))**2*aIpha**2*yjl*yjl/pi/np 


do  170  k=l,ns 

scahh(k)=c*scahh(k) 

scavh(k)=c*scavh(k) 

scahv(k)=c*scahv(k) 

170  scavv(k)=c*scavv(k) 

c— derivatives  of  qext,cvv  and  chv/vertical  inddence/ 
esiim=(0.0,0.0) 

fsuni=(0.0,0.0) 

do250j=l^p 

cin=anplx(cos(z(j)),-sin(z(j))) 

esuin=esuni+  (e(l  j)*2.0*ifindx+{epslon-1.0)*el(l  j))*cni 
fsum=fsuin+  (ex(lj)-cinplx(0.0^(j))*e(lj)/alpha)*cm 

250  continue 

qetn=4.0/np*yj  1  *ainiag(esum) 
qetk=4.0/np*yj  1  *real(esnm) 

qetx=(1.0^.3*alpha2)/alpha/(1.0-0.1*alpha2)*qext 
&  +  4.*alpha/3.0/np*aimag((epslon-1.0)*fsum) 

temp0=(epslon-1.0)*conjg((epslon- 1 .0)) 
tenipl=alpha2*yj  1  *yj  1/pi/np 

tempn=4.0*real(rfindx)*(rfindx*conjg(ifmdx>l  .0) 
tempk=4.0*ainiag(rfindx)*(riindx*conjg(rfindx)+l-0) 

do230k=Ms 
cangIe=cosang(k) 
sangle=sinang(k) 
esuni=(0.0,0.0) 
fsuin=(0.0,0.0) 
ep(l)=0.0 
ep(2)=0.0 
cp(5)=(00.0.0) 
cph(5)=(0.0,0.0) 
do220j=l,np 
ij=z(j)*cangle+y0)*sangle 
cm=cmpIx(cos(ij),-sin(ij)) 
esum=esuni+  e(l  j)*an 
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fsuin=fsum+  el(l  j)*cm 
ep(l)=ep(l)+(cangle*e(2j)-sangle*e(3j))*cni 

cp(2)=ep(2)+{caiigle*el(2j)-sangle*cl(3j))*cm 

ep(5)=ep(5>Kex(l  j)-d*ij/alpha*e(U))*cm 

eph(5)=eph(5)+(cangIc*(ex(2J)-c(2j)*ci*i3/alplia) 

&  -sangle*(ex(3j)-e(3j)*ci*ii/alpha))*cin 

220  continue 

flicta=esum*conjg(fsuni)4conjg(esuin)*fsuni 

qj(3)=ep(l)*conjg(ep(2)>HXHijg(ep(l))*cp(2) 

phi=cmplx(0.0,1.0)*(conjg(esuin)*fsuni-esuin*conjg(fsum)) 

cp(4)=anplx(0.0,1.0)*(conjg(ep(l))*ep(2><p(l)*conjg(ep(2))) 

dkvvn(k>=scaw(k)/tenipO*tempn+tempO*templ  *theta 

dkwk0t)=scaw(k)/tempO*tempk+tempO*tenipl*phi 

dkhvn(k)=scahv(k)/tenip0*tempn+temp0*templ*ep(3) 

dlchvk(k>=scahv(k)/tempO*tempk+tempO*templ  *ep(4) 

dkvvx(k)=(4.0-1.2*alpha2yalpha/(l-0.2*alpha2)*scaw(k) 

&  +temp0*templ*(ep(5)*conjg(esum>+conjg(ep(5))*esuni) 
dkhvx(k)=(4.{)-1.2*alpha2)/alpha/(l-0.2*alpha2)*scahv(k) 

&  -rteinp0*templ*(eph(5)*conjg(ep{l)>«onjg(eph(5))*ep(l)) 
230  continue 


c — derivatives  of  qexth,  cvh,  &  chh 
esuinh=(0.0,0.0) 
fsunih=(0.0,0.0) 

do350j=l4ip 

cm=cmplx(cos(zO')),-sin(z(j))) 

csunih=csuinh+(ch(2j)*2.0*ifindx+(epslon-1.0)*ehl(2J))*cm 

fsuinh=fsunih+(ehx(lJ)-anpIx(O.OXj))*eh(lj)/alpha)*cni 

350  continue 

qehtn=4.0/np*yj  1  ♦aimag(esumh) 
qehtk=4.0/np*yj  1  *real(esuinh) 

qehtx=(l.(W)-3*alpha2)/alpha/(l,(M).l*alpha2)*qexth 
&  +  4.*alpha/3.0/np*ainiag((epslon-1.0)*fsuinh) 

temp0=(cpslon-1.0)*conjg((epslon-1.0)) 
teinpl=alpha**2*yjl*yjl/pi/np 
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temim=4.0*real(ifindx)*(ifindx*conjg(ifindx)-l  .0) 
tcmpk=4.0*aimag(rfiiidx)*(ifindx*conjg(ifindx)+-1.0) 
do310k=l,ns 
cangle=cosang(k) 
sang^e=aiiang(k) 
esuinli=(0.0,0-0) 
fsuiiih=(0.0A0) 
epKl)=0.0 
eph(2)=0.0 
ep(5)=(0A0.0) 
cph(5)=(0A0.0) 
do320j=l^p 
ij=z(i)*cangle+y(j)*sanglc 
cin=cmplx(cos(ij),-sin(ij)) 
esuinh=esumh+  eh(I  j)*cm 
fsuinh=fsumh+  ehl(l  j)*cin 

ei*(l)=cph(l>+<cangle*ch(2j)-sangle*eh(3j))*cm 

eph(2)=eph(2>f(cangle*ehl(2j)-sangle*ehl(3j))*cm 

q)(5)=ep(5>Kel«(lj)-«*d/alpha*eh(lj))»cni 

e|A(5)=eph(5>t<cangle*(elix(2j)-ch(2j)*ci*ij/alpha) 

&  -sangle*(e!ix(3j)-ch(3j)*ci*ij/alpha))*cm 

320  continue 

tiietah=esumh*conjg(fsumh)+conjg(csunih)*fsumh 

eph(3)=eph(l)*conjg(ci*(2)>+conjg(cph(l))*eph(2) 

pluh=cniplx(0.0,1.0)*(conjg(esunih)*fsuinh-esuinh*conjg(fsumh)) 

eidi(4)=cmplx(0.0,l-0) 

$  *(conjg(cph(l))*cph(2>eph(l)*conjg(eph(2))) 

dkvhn(k)=scavh(kytempO*tempn+tempO*tenipl*thetah 
dkvhk(k)=scavh(k)/tempO*tempk+tcnipO*templ*phih 
dkhhii(k)=scahhfkVtempO*tempn+tenipO*templ  *eph(3) 
dklihk(k)=scahh(k)/tenipO*tempk+tenipO*tenipl  *eph(4) 
dkvhx(k)=(4.0-1.2*alplia2)/alpha/(l-0.2*alpha2)*scavh(k) 

&  +temp0*templ*(ep(5)*conjg(esuiiih)+conjg(ep(5))*esunih) 

rfifhhT(k>=(4  n-1 ,2*alpha2Valpha/(  l-0.2*alpha2)*scahh(k) 

&  +tcmpO*templ*(eph(5)*conjg(eph(l)>+conjg(eph(5))*eph(I)) 
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310  continue 


c — orientation  averages 

teinp=w(ic)*w(ip)*chisin*pi/8.0 

qexl0=qext04qext*temp 

qabsO=qabsOfqabs*ten]p 

qscaOaqscaOfqsca*ten^ 

qexthO=qext]i&fqextfa*temp 

qabshO=qabshOtqabsh*temp 

qscali0=qscah04qscah*teinp 

qetiiO=qetnO+qetn*temp 

qetkO=qetkOfqctk*temp 

qehtnOaqehtnOfqehtn^temp 

qehtlc0=qehtkO4qehtk*temp 

qetxO=qctxOfqetx*teinp 

qehtxO=qehtxO+qeht!c*temp 

do360k=l,ns 

scavvO0c)=scavvO(k)+scaw(k)*temp 

scahvO(k)=scahvO(k)+scahv(k)*temp 

scahhO(k)=scahhO(k)+scahh(k)*temp 

scavhO(k)=scavhO(k)+«cavh(k)*tcmp 

dkvvnO(k)=dkvvnO(k)+dkvvn(k)*temp 

dliwkO(k)=dkvvkO(k)+<Ikvvk(k)*teinp 

dkhvnO(k)=dkhvnO(k>tdkhvn(k)*teinp 

dkfavkO(k)=d1divkO(k>Klkhvk(k)*temp 

dklilm0(k>=dkfaim0(k)4<lkhlm(k)*tem 

(]khhkO(k)=dkhhkO(k>t<Ikhhk(k)*ten^ 

<lkvlinO(k)=dkvhnO(k)tdkvhn(k)*temp 

dkvhkO(k)=dkvhkO(k)^<lkvhk(k)♦te^^ 

dkvvxO(k)=dkvvxO(k)+dkvvx0r)*temp 

dkhvxO(k)=dkhvxO(k)tdkhvx0t)*temp 

dkvhxO(k)==dkvhxO(k>t<lkvhx(k)*temp 

dkhhxO(k)=dkhhxO(k)4dkhhx0t)*teinp 


360  continue 


410  continue 
420  continue 


return 

end 


c¬ 


c — agglouO.  orientation  averaged  version  of  aggloul 

c - gaussian  integration  (ng  points,  ng  <=20)/3-12-95 - 

c - — - ^wlou  - 

subroutine  agglouO(np,alpha,epslon,ns) 
parameter  (innp=l(X)) 

complex  epslon,e(3,mnp),esum/sum,ep(6),theta,phi, 

$  cm,  eh(3,mnp),  eph(6),  esumh4sumh,thetah,phih,ci, 

$  el(34imp),ehl(3,mnp),ex(3,mnp),ehx(3,mnp),  ctem 
character  pol*  1  .chain*  1  ,cff*  1 
real  scahh(181),scavv(181),scahv(181),scavh(181) 
common/p2/pol,efr,chain 
common/lou/lu,lu  1  ,lu2,Iu3,iu4 

common/loul/cospsi(20),sinpsi(20),coschi(20),sinchi(20),w(20),ng 

common/jl/e,  di 
common/jI2/el,  ehl 
common/jI3/ex,ehx 
common/12/x(mnp),y(mnp),z(mnp) 
common/1  l/xO(mnp),yO(mnp),zO(mnp) 
c — orientation  aveiagedquantities 

conunon/j2s/scahh0(181),  scavv0(181),  scahv0(181),  scavh0(181), 
$  cosang(181),  sinang(181) 

common/qqcxt/qextO,qabsO,qscaO,qexthO,qabshO,qscahO 

pi  =3.1415926 

pi2=pi*0.5 

alpha2=a]pha*a]pha 

ci=cmplx(0.0,1.0) 

ifGu.eq.2)then 

yjl=(sin(alpha)/aIpha-cos(alpha))/alpha 

else 

yjl=alpha/3.0 
end  if 

gang=4.*yjl/np 

c=(epslon-l.)*conjg((epslon-1.0))*alpha2*yjl*yjl/pi/np 
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call  coord(np^pha,l.A04*A0) 
call  zhoul(np,alpha,epslon) 

c— data  initialization - 

qext0=0. 
qsca0=0. 
qabs0=0. 
qcxth0=0. 
qscah0=0. 
qabsh0=0. 
do  10  k=:l,ns 
scahh0(k)=0. 
scahv0(k>=0. 
scavv0(k)=0. 
scavh0(k)=0. 

10  continue 
do  420  ip=Mg 
psicos=cospsi(ip) 
psisin=sinpsi(ip) 
do410 

chisin=sinchi(ic) 

chicos=coschi(ic) 

call  coord(np,alpha,chicos,chisin,psicos,psisin) 

if(eff,cq.'n*  .or.  eff.eq.V)then 
call  zhou0(np,alpha,cpslon) 
else 

call  2hou3(np,alpha,epslon) 
end  if 

c — extinction  &  absorption  cffiency  factors 
asun3=0.0 
csuni=(0.,0.) 
asumh^O.O 

esumh=(0.,0.) 

do40j=Mp 


ctem=cmplx(cos(z(j)),-sin(z(j))) 
esum=esum+ e(I  j)*ctem 

asuin=asuin+e(lj)*conjg(e(l  j)>+€(2J)*coijjg(e(2j)) 

$  4«(3o)*conjg(c(3j)) 

esiiinh=esunih+eh(2j)*ctem 

asunih=asuinh+eh(lj)*conjg(eh(lj)>+€h(2j)*conjg(eh(2j)) 

$  4€h(3J)*conjg(eh(3j)) 

40  continue 

qext=aiinag((epslon-l.)*esuin) 
qabs=ainiag(epsIon- 1 .0)*asum 
qexth=ainiag((epslon-l.)*esunih) 
qabsh=aiinag(cpslon- 1 .0)*asuinh 
qext=gang*qext 
qabs=gang*qabs 
qexth=gang*qexth 
qabsh=  gang*qabsh 

c— total  scattering  effiency  factor,  set  lul=0  to  skip 
ifGul.eq.O)then 
qsca=qext-qabs 
qscah=qexth-qabsh 
else 

call  louO(np, alpha, epsion,qsca,qscah) 
end  if 

do  140  k=l,ns 
cangle  =  cosang(k) 
sangle  =  sinang(k) 

dol20I=U 

ep(l)=cmpix(0.,0.) 

120  ephO)=cmplx(0.,0.) 

dol32j=Mp 

ij=z(j)*cangle+y0*sanglc 
ctein=cmplx(cos(ii),-sin(ij)) 
do  130 1=U 
epO>=ep(l)+cteni*e(lj) 
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130  cphO)=eph(l>fctem*eh(lj) 

132  continue 

theta=cangle*ep(2)-sangle*ep(3) 

plii=cp(l) 

scahv(k)=theta*conjg(thela) 

scaw(k>=phi*conjg(phi) 

thetah=cangle*eph(2)-sangle*eph(3) 

phili=eph(l) 

scahh(k)=thetah*conjg(tiietah) 

scavh(k)=phih*conjg(phih) 

140  continue 
do  170  k=l,ns 

scalih(k)=c*scahh(k) 

scavh(k)=c*scavh(k) 

scahv(k)=c*scabv(k) 

170  scaw(k)=c*scavv(k) 


c — orientation  averages 

temp=w(!c)*w(ip)*chisin*pi/8.0 

qextO=qextOfqext*temp 

qabsO=(]absO+qabs*temp 

qscaO=qscaOfqsca*lenip 

qexthO=qexthO+qexth*temp 

qabshO=qabshOfqabsh*tenip 

qscah0=qscah04qscab*tenip 

do  360  k=Ms 

scavv0(k>=scavv0(k>«caw(k)*temp 

scahv0(k)=scahv0(k>+scaliv(k)*temp 

scal»h0(k)=scahh0(k)+scahh(k)*tenip 

scavh0(k)=scavh0(k)+«ca>1i(k)*temp 


360  continue 
410  continue 
420  continue 


return 

end 
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c - - - - - - 

c— sub:  zhouO  solves  the  internal  field  under  non-coupling  option 
c —  (efl='n')  in  the  measuring  coordinate  system  (x,y,z) 

c~  it  also  gives  the  internal  field  of  the  raylcigh-gans  theory 
c_  (eff=V)  adikdonfiiday,  10-14-94  by  w.  lou 

c - - — - - — - — 

subroutine  zhouO(np,alpha,epslon) 
parameter  (mnp=:  100,  mn3=mnp*3) 
complex  epslon,b(mn3),bh(mn3),c,d,cl 
character pol*l,  efT*l,  diain*l 
commonyp2/pol,eff,chain 
commonyiou/lu,Iu  I,lu24u3,lu4 
common/jl/b.bh 

common/12/x(mnp),y(mnp),z(mnp) 
if(eff  .eq.  'r')then 
do  60  i=l,np 

c=cmplx(cos(z(i)),sin(z(i))) 

niowO=3*C>-l) 

bh(nrowOf  1  MO.0,0.0) 

bh(nrowOt-2)=c 

bh(nrowOf3M0.0,0.0) 

b(nrowO<-l)=c 

b(nrow0+2)=(0.0,0.0) 

b(nrow0+3)=(0.0,0.0) 

60  continue 
else 

cl=3./(epslon+2.) 

if0u.eq.l)then 

dl=2.0*alpha/3.0 

d=(1.0,0.>-alpha*alpha*(epslon-l.)/(epslon+2.)*cmplx(1.0,dl) 
else  if(lu.eq.0)then 
d=(1.0,0.0) 
else  if(lu.eq.2)then 

yjl=(sin(alpha)/alpha-cos(alpha)yalpha 

yyl=-{cos(alpha)/alpha-win(alpha))/alpha 
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yifa  1 =cinplx(-yy  1  ,yj  1 ) 

d=(3.0*epslon-(e‘psIon-l.)*2.*alpha**2*yihl)/(epslon+2.) 
else 
end  if 

do64i=l4ip 

c=cmplx(cos(z(i)),sin(z(i))) 

nrowO=3*(i-l) 

bh(nrow0+l)=(0.0,0.0) 

Wi(nrow0+2)=c*cl/d 

bh(nrow0f3)=(0.0,0.0) 

b(nrowOfl)=c*cl/d 

b(nrow(>+2)=(0.0,0.0) 

b(nrow0f3)=(0.0,0.0) 

64  continue 
end  if 
return 
end 
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c— sub:  zhoul  establishes  the  matrix  a  for  the  internal  field 
c —  in  the  structural  coordinate  system  (x’,y'^') 
subroutine  zhoul(np,alpha,epslon) 

parameter  (mnp=  100,  nin3=mnp*3) 

complex  epslon,a(rtur3,rrm3),  al(mn3,mn3),  t(33)jC*d,ctem 

character  efI*l,pol*l , chain*  1 

commorr/lou/lu,lu  1  ,lu2,lu3  ,lu4 

common/jla/a,al 

common/1  l/x0(rrmp),y0(mnp),z0(mnp) 
common/p2/poI,e£F,chain 
pi=3. 1415926 
nrank=3*np 
if  (cff.eq.'n')then 
return 

else  if  (eff-eq-f )  then 
ctem=cmplx(0.,-l.) 
else  if  (eff.eq.'w')  then 
ctem=cmplx(0.,l.) 
oidif 

c=alpha**3*(epslon-l.)/(3.*(epslon+2.))*ctem 

ifOu.cq.l)then 

dl=2.0*alpha/3.0 

d=(10,0.)-alpha*alpha*(epslon-l.)/(epslon+2.)*cmplx(1.0,dl) 
else  if(lu.eq.O)then 
d=(l  .0,0.0) 
else  rf(lu.eq.2)then 

yj  l=(sin(alpha)/alpha-cos(alpha))/alpha 
c=alpha**2*(epslon-l.)/(epslon+2.)*ctem*yjl 
yyl=-(cos(alpha)/alpha+sin(alpha))/alpha 
yihl=cmplx(-yyl,yjl) 

d=(3.0*epsIon-(epslon-l.)*2.*alpha**2*yihl)/(epslon+2.) 

else 
end  if 

dol0i=l4ip 
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nrow(>=l+3*0-l) 

im)we=3*i 

do  20  iirow=iiix)wO,nrowe 
do  20  ncol^irow^we 
a(niDw^col)=(0-0.0-0) 

if  (ncol.eq.iirow)then 
a(iirow^col)=d 

adif 

20  continue 

if(i  .eq.  np)  goto  10 
ir=i+l 

do30  j=ir4)p 
dx=x0(j)-x0(i) 
dy=y0(j>y0Ci) 
dw00>z0(i) 

dj=sqit(dx*dx+dy*dy+dz*dz) 

dj2=sqrt(dx*dx-Hly*dy) 

ducos=dz/dj 

chisin=dj2/dj 

if(dj2.1e.0.1c-30)then 

psicos=1.0 

psisin=0.0 

else 

psicos=dx/dj2 
psisin=dy/dj2 
end  if 

djk=dj*2.0*alpha 

call  tmatrx(djk,chicos,chisin,psicos,psisin,t) 

ncol0=3*a-l) 

do  40  nc=l,3 
do  40  nrow=jm)w04irowc 
40  a(nrowjncol0fnc)=c*t(nrow-nrow0fl4ic) 
30  continue 
10  continue 


do  50  niow=2,niank 


iicddg=nrow-l 
do  50  ncoI=l^coldg 
50  a(iirow4icol)=a(ncol4in)w) 
if(eff.cq.'f)then 
do  56  i=l,iirank 
do56  j=l,iiraiik 
56  al0d)=a(ij) 

call  cbdecp(a,nrank) 
else 

do  60  i=l4iTaiik 
60  a(i,i)=(1.0,0.0) 

emfif 


retUTD 

end 


c- 


c— sub:  zhou3  solves  the  internal  field  in  the  (x',y',z')  system 
c —  then  transforms  it  into  the  measuring  coordinate 

subroutine  zhou3(np,alpha,epslon) 
parameter  (innp=l(X),  mn3=mnp*3) 
complex  epslon,a(™D34nn3),b(inn3),bh(iiui3),c,d  ,bl(mn3), 
$  ex,  ey,  ez,  exh,  eyh,  czh,  al(mn3,mn3), 

$  bs(inn3),  bhs(mn3) 
realt(33) 

diaiacter  efPl,pol*l,chain*l 
common/jl/b,bh 
conunon/jle/bs,bhs 
common/jla/a,al 

common/12/x(mnp),y(mnp)Xnu>P) 

common/p2/pol,efr, chain 

common/lou/lu,lul,lu2,lu3,iu4 

common/transl/t 

nrank=3*np 

do64i=l4)p 

c=cmplx(cos(z(i)),sin(z(i))) 

niDwO=3*(i*i) 

bh(nrow0r-l)=c*t(2,l) 

Wi(im)w(>f2)=c*t(2,2) 

bh(nrowO+3)=c*t(23) 
b(niowOr-l)=c*t(l,l) 
b(nit)w(>f2)=c*t(l,2) 
b(niDwO+3)=c*t(l  3) 

64  continue 

d=37(epslon+2.) 
if(eff.eq.f)then 
do70  i=l,iirank 
bC«)=bCi)*d 
70  bh(i)=bha)*d 

call  cbsolx(a,b,nrank) 
call  cbsolx(a,bh,nrank) 
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else  if(eflf.eq.’w')tbcn 
do  80  i=l4iraDk  • 
bl(i)=cmplx(0.,0.) 

do80  j=l,i)rank 
80  bl(i)=blCiHaa4)*bO) 
do  85  i=:l,iiraiik 
b(i)=cmplx(0.,0.) 

do85  j^l^uank 
85  Kj)=b(i)+a(j4)*bh0) 

if(Iu.eq.l)then 
dl=2.0*aIpha/3.0 

c=(1.0,0.)-alpha*alpha*(epslon-l.)/(epslon+2.)*cmplx(1.0,dl) 
else  if0u.eq.0)theii 
c=(l. 0,0.0) 
else  if(lu.eq.2)then 

yj  1  =(sin(alpha)/aIpha-cos(alpha))/aIpha 

yyl=-(cos(alpha)/alpha+sin(alpha))/alpha 

yihl=cmplx(-yyl,yjl) 

c=(3.0*epslon-(epslon-l.)*2.*a]pha**2*yihI)/(epsloD+2.) 

else 
end  if 

do  90  i=l4irank 
bh(i)=b(i)*d/c 
90  b(i)=bl(i)*d/c 
else 
stop 
end  if 

do  100  i=l,np 

il=3*(i-lKl 

12  =  il  +  l 

13  =  il  +2 
bs(il)=bai) 
bsCi2)=bCi2) 
bs(i3)=b(i3) 
blis(il)=  Wi(il) 
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bhs(i2)=  bh(i2) 
bhs03)=bh(i3) 

ex=b(il)*t(l,l)+b(i2)*t(U)+b(i3)*t(U) 
ey=b(il)*t(2,l)+b(i2)*t(2^>fb(i3)*t(23) 
cz=b(il)*t(3,l)+Ki2)*t(3^)+ba3)*t(33) 
e3di=bh(il)*t(l,l)+bh(i2)*t(U)+bha3)*t(U) 
cyh=bh(i  1  )*t(2,l  )+bh(i2)*t(2^)+bh(i3)*t(2,3) 
ezh=bh(il)*t(3.1>+bh(i2)*t(3^>+bh(i3rt(33) 
b(il)=ex 
ba2)=ey 
bCi3>=ez 
bh(il)=exh 
bh02>=eyh 
UiQ3)=ezh 
100  continue 
return 
end 
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c- 


c— sub:  louO  computes  total  scattering  efficiency  factor 
c —  in  the  structural  coordinate  system 

subroutine  louO(np,alpha,epslon,qsca,qscah) 
parameter  (mnp^lOO,  nm3=mnp*3) 
complex  epslon,a(™t3,mn3),c,d  ,al(mn3,mn3), 

$  bs(nm3),  bhs(mn3),bji0,bjlh0,ctem 
character  cff*'  1  ,pol  *  1  ,chain*  1 
common/p2/pol,eff,chain 
common/jle/bs,bhs 
common/jla/a,al 
common/lou/lu4ul,Iu2,lu3,lu4 
if  (eff.eq.*n7then 
return 

else  if  (eff.eq.T)then 
ctem=cmplx(0.,- 1 .) 
else  if  (eff-eq-V*)  then 
ctem=cmplx(0.4.) 
endif 

c=alpha**3*(cpslon-l.)/(3.*(epslon+2.))*ctem 

yjl  =  alpha/3.0 
ifOu.cq.l)then 
dl=2.0*alpha/3.0 

d=(10,0.)-alpha*alpha*(epslon~l.)/(cpslon+2.)*cmplx(1.0,dl) 
else  ifOu.eq.O)then 
d=<L0,0.0) 
else  if(Iu.cq.2)then 

yj  l=(sin(alpha)/alpha-cos(alpha))/aIpha 
c=aIpha**2*(epsIon-l.)/(epslon+2.)*ctem*yjl 
else 
endif 
bjl=0.0 
bjlh  =  0.0 
do64j=l,np 
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do  64 1=1, np 
if(j  .eq.  l)then  ' 
do50j2  =1,3 
jl  =(j-l)»3+j2 

bjl  =bjl  +  2.0*  bs(il)*conjg(bs(jl)) 
bjlh  =bjlh  +  2.0*  bhs(jl)*conjg(bhs(jl)) 
50  continue 

else 

do  54  j2  =1,3 

jl=0-l)*3+j2 

bjlO  =(0.0,0.0) 

bjlhO  =(0.0,0.0) 

do  52 12  =13 

11=(1-1)*3+12 

z  =  real(aiai41)/c) 

bjlO  =  bjlO  +  z*conjg(bs(ll)) 

52  bjlhO  =  bjlhO  +  z*conjg(bhs01)) 

bjl  =bjl  +bjlO*bsai) 

54  bjlh  =bjlh  +bjlhO*bhs(jl) 
end  if 

64  continue 

qsca  =  bjl*4.0/3.0/np*alpha*alpha*yjl*yjl 
$  *(epslon-1.0)*conjg(epslon-1.0) 
qscah  =  bjlh/bjl*qsca 
return 
end 
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c¬ 


c — sub:loul 

c —  generate  cosine  and  sine  of  the  gausian  points(ng  points 
c - dateof  generation:  december21,  1994  - lou- 


subroutine  ioulO 
parameter  (ngs=10) 

rBaIzx(20),w(20)^6(3),w6(3)^8(4),w8(4);KlO(5),wl0(5), 

&  zxI2(6),wl2(6)^16(8),wl6(8),zx20(10),w20(10) 

coiiimon/loul/cospsi(20),sinpsi(20),coschi(20),sinchi(20),w4ig 

data  pi/3.1415926/ 

data  2x16/0.09501251,0.28160355,0.45801678,0.61787624, 

&  .75540441,.86563120,.94457502,.98940093/,  wl6/0.18945061, 
&  0.18260342,0.16915652,0.14959599,0.12462897,0.09515851, 
&  0.06225352,0.02715246/ 

data  zxl2/.12523341,.36783150,.58731795,.76990267,.9041 1726, 
&  .98156063/,  wl2/.24914705,.23349254,.20316743,.16007833, 
&  .10693933,.04717534/ 

datazxl0/.14887434,.43339539,.67940957,.86506337,.97390653/, 
$  wl0/.29552422,.2692667,.2190864..14945135,.06667134/ 
data  2x8/  .18343464,.52553241,.79666648,.96028986/, 

&  w8/.36268378,.31370665,.22238103,.10122854/ 
data2x6/.23861919,.66120939,.93246951/, 

&  w6/  .46791393,0.36076157,0.17132449/ 
pi2=pi*0.5 
ng=ngs 
if(ngs  .eq.  6)tben 
do  12  i=l,ngs/2 

2x(i+ngs/2)=^6(i) 

2x(i)=  -zx6(ngs/2-i+l) 
w(i+ngs/2>=w6(i) 

12  w(i)=w6(ngs/2-i+l) 

else  if(ngs  .eq.  8)then 
do  13  i=l4igs/2 

2x(i+ngs/2)=2x8(i) 
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do20  i=l,ng 
psi=(zx(i>i-l.)*pi 

cospsi(i)=cosQ)si) 
smpsiCi)=sin(p8i) 
20  continue 


do40i=l4>g 
chi=(zx(i)+l-)*pi2 
sinchi0)=sin(chi) 
coschi(i)=cos(chi) 
40  continue 
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return 


end 
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c 


c — sub:zhou31 

c - zhou3  with  derivatives  of  intenial  field  wrt  (n4c)&  xp 

c -  output  in  conunon/jl2/&/jl3/.  full  coupling  only. 

c - ^lou,  august  20(n,k),  and  September  28,  1994  for  xp's  - 

subroutine  zhou31(np,alpluMfindx) 
parameter  (innp=100,  mn3=mnp*3) 

complex  rfindx,epsIon,a(mn3,mn3),b(mn3),bh(mn3),c,d,cl,clx, 

&  ex,  ey,  ez,  cxh,  eyh,  ezh,  bl(mn3),bhl(mn3), 

&  al(mn3,mn3),  bx(mn3),bhx(mn3),h0,hl,li2,bs(mn3),bhs(inn3) 

realt(33) 

character  eff  1  .chain*  1  ,pol  *  1 

common/jl/b,bh 

common/jle/bs,bhs 

common/jI2/bl,bhl 

common/jl3/bx,bhx 

common/jIa/a,al 

common/12/x(mnp),y(mnp),z(mnp) 
common/p2/pol,efr, chain 
common/lou/lu,lul,lu2,lu3,lu4 
common/transl/t 

nrank=3*np 

epsl<ai=rfindx*ifindx 

ez=(®pslon-1.0)/(epslon+2.0) 

d=3./(epslon+2.) 

a]pha2=a]pha*alpha 

sO=sin(alpha) 

cO=cos(aIpi)a) 

zl=:c(M-alpha*sO 

z2=s0-a]pba*c0 

cl=d4cz*(3.0-2.0*cmplx(zl,z2)). 

clx=-2.0*cmplx(c0,s0)*ez 

c — incoming  wave  rotation - 

do64i=l,np 
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c=cmplx(cos(z(i)),sin(z(i))) 

nrowO=3*G-l)  ' 

bh(mow0+l)=c*t(2,l) 

bh(mow(H-2)=c*t(2^) 

bh(iirowO+3)=c*t(23) 

b(im)w04-l)=c*t(l,l) 

b(iirow(>f2)=c*t(l 

b(im)wO+3)=c*t(l  3) 

64  continue 

c — internal  field  in  the  structural  cowdinate  system,  including 

c - derivatives  of  the  field,  full  coupling  only. 

if  (efT.eq.T)then 
do  70  i=l,nrank 
b(i)=bCi)*d 
bh(i)=bh(i)*d 

bia)=bC0 

70  bhl(i)=bh(i) 

call  cbsolx(a,b4irank) 
call  cbsolx(a,bh,nrank) 

c - derivatives  w.r.t  size  parameter  alpha 

c=cmplx(0.0,-1.0)*alpha2*alpha*ez/3.0 
do  58  l=l,np 
do56m=l,3 
nl=3*0-I)+ni 

ex=(0.0,0.0) 

cy=(0.0,0.0) 

exh=(0.0,0.0) 

eyh=(0.0,0.0) 

do54j=l,np 

n2=3*0-l)+m 
n3=3*(j-l) 
if(j.eq.l)go  to  54 

iho=sqit((xa)-x(j))*(xa>xO))+<yO>-y(j)) 

$  *(yO)-y(j)Mz(l>z(j))*(z(l)-za))) 

s0=sin(iho) 
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cO=cos(tho) 

hO=cmpIx(sO,-cO)/ito 

zl=(3./rho/riio-l.)*sO-3.*cO/rho 

z2=(l.-3./riio/rho)*cO-3.*sO/iiio 

|]2=cmplx(zl  ^yibo 

hl=rho/3.0*(h(>fli2) 

ex=ex  +  2.*c*b(n2)*(rho*hl  -  3.*h0  +  rho*hl*h0/h2) 
ey=ey+  iho*hl/h2*(al(nl^+l)*b(n3+l) 

&  +  al(nU3+2)*Kn3+2Hal(nl,ii3+3)*b(ii3+3)) 

exh=exh  +  2.*c*bh(n2)*(rbo*hl  -  3.*h0  +  ilio*hl*h0/li2) 
eyh=eyh  +  rtio*hl/h2*(al(nl^3+l)*bh(n3+l) 

&  +  al(nM3+2)*bh(n3+2)+al(nl^+3)*bh(n3+3)) 

54  continue 

bx(nl)=bl(nl)*cmplx(-alpha/3.0^yalplia) 

&  +  b(nl)*alpha*(cl/3.0<lx)+  (ex-eyyalpha 

56  blix(nl)=bhl(nl)*cmplx(-alpha/3.0,z(iyalpha) 

&  +  bh(nl)*alpha*(cl/3.0-clx)+  (exh-eyhyalpha 

58  continue 

call  cbsolx(a,bx4irank) 
call  cbsolx(a,bhx4iiank) 

c - derivativesw.r.t  refractive  index  (n,k) 

do  72  i=l4uank 

bhl(iMd*bh(i)-bhl(i)y(epslon-1.0)*2*ifmdx 
bl(i)=(d*bCi)-bl(i)y(epslon-1.0)*2.0*ifindx 
72  continue 

call  cbsolx(a,bl4uank) 
call  cb8olx(a,bhl,nrank) 
else  if(eff.eq.’w')then 
do80i=l4uank 

bl(i)=cmplx(0.,0.) 

do  80  j=l4iiank 
80  blCi)=blCiHa(jJ)*b(j) 
do  85  i=l4uank 
b0)=cmplx(0.,0.) 
do85  j=l4irank 
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85  b(i)=b(i)+aa4)*bh(j) 

if(Iu.eq.l)then 
dl=2.0*aIpha/3.0 

c=(lA0.>alpha*aIpha*(epslon-l.)/(epsIon+2.)*cmplx(1.0,dl) 
else  ifOu.eq.O)then 
c=(1.0,0.0) 
else  if(lu.eqJ2)then 

yjl=(sin(alpha)/alpha-cos(alpha))/alplia 

yyl=-{cos(alpha)/alpha+sin(alpha))/alpha 

yihl=cmplx(-yyl,yjl) 

c=(3.0*epslon-{epslon-l.)*2.*alpha**2*yihl)/(epslon+2.) 
else 
end  if 

do  90  i=l,nrank 
Mia)=b(i)*d/c 
90  b(i)=bl(i)*d/c 
else 
stop 
end  if 

c—field&derivative  rotate  back  to  the  measuring  coordinate  system 
do  100  i=l4»P 
il=3*(i-l)+l 
i2  =  il  +  1 
i3=il  +  2 
bs(il)=b(il) 
bsa2)=ba2) 
bs(i3)=b(i3) 
bhs(il)=bh(il) 
bhs(i2)=  bh{i2) 
bhs(i3)=bh(i3) 

ex=b(il)*t(l,l)+b(i2)*t(U>4-b(i3)*t(U) 

ey=b(il)*t(2,l)+b(i2)*t(2;i>fb(i3)*t(23) 

ez=bai)*t(3,lHba2)*t(3^>+Ki3)*t(3,3) 

exh=bhai)*t(l.l)+bh(i2)*t(UKbha3)*t(13) 

eyh=bhai)*t(2,lHbh(i2)*t(2^>+bh(i3)*t(23) 
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ezh=bh(il)*l(3,l)+bh(i2)*t(3^>fbh(i3)*t(33) 
b01)=ex 
bCi2)=ey 
bfi3)=ez 
U)01)=%xh 
bha2)=cyh 
l*(i3)=ezh 
100  continue 
do  102  i=l^p 

11  =3*(i-l)+l 

12  =  il  +  1 
i3=il  +  2 

ex=bl(il)*t(l,l)+bl(i2)»t(U)+bl(i3)*t(13) 
ey=bl01)*t(2,l)+bia2)*t(23)+bl(i3)*t(23) 
cz=bl(il)*t(3,l)+bl(i2)*t(33)+bl(i3)*t(33) 
cxh=bhl(il)*t(l,lHbhl(i2)*t(13)+bhia3)*t(13) 
cyh=bhl(il)*t(2,l)+bhl(i2)*t(23)+blil(i3)*t(2,3) 
czh=bhl(il)*t(3,l)+bhia2)*t(33)+bhl(i3)*t(33) 
blCil)=ex 
blCi2)=ey 
bl03)=ez 
bhlfil)=exh 
bhl(i2)=eyh 
bhlfi3)=czh 
102  continue 
do  104  i=l^p 

11  =3*(i-l)+l 

12  =  il  +  1 
i3=il  +  2 

ex=bx(il)*t(l,l)+bx(i2)*t(U)+bx(i3)*t(l,3) 

ey=bxai)*t(2,l}+bx(i2)*t(23)+bxa3)*t(23) 

ez=bx(il)*t(3,l)+bx(i2)*t(33)+bx(i3)*t(33) 

exh=bhx(il)*t(l,l)+bhx(i2)*t(13)+bhx(i3)*t(13) 

eyh=bhx(il)*t(2,l)+bhx(i2)*t(23)+bhx(i3)*t(2,3) 

ezh=bhx(il)*t(3,l)+bhx(i2)*t(33HblK03)*t(33) 
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bx(il)=cx 
bx02)=ey 
bx(i3)=ez 
bhx01)=exh 
bhx(i2)=eyh 
bhx03)=ezh 
104  continue 
return 
end 
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c 


c — subrtmatrx 

c - calculates  matrix  t 

c 

subroutine  tinatrx(z,cc,sc,cp,sp,t) 
complex  t(33)>bl0,hl2 
s=^(z) 

(?=cos(z) 

hlO=cmplx(s,-cyz 

zl=(3./z/z-l.)*s-3.*c/z 

z2=(1.-37z/z)*c-3.*s/z 

hl2=cmplx(zl;z2)/z 

t(l,l)=2.0*hl0-hl2*(l-3.0*sc*sc*cp*cp) 

t(l,2)=hl2*3.0*sc*sc*sp*cp 

t(U)=hl2*3.0*sc*cc*cp 

t(2,l)=t(U) 

t(2;i)=2.0*hl0-hl2*(1.0-3.0*sc*sc*sp*sp) 

t(2,3)=hl2*3.0*sc*cc*sp 

t(3,l)=t(U) 

t(3^)=t(2,3) 

t(33)=2.0*hl0-hl2*(l-3.0*cc*cc) 

return 

end 
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c 


c — sub:  strait 

Q - generate  straight  chain  geometry 

subroutine  strait(nt) 

parameter  (mnp=100) 

common/1  l/x(ninpXy(mnp)Xiimp) 

common/ran/xl 

nt2=nt/2 

if(nt-2*nt2.eq.0)then 
do  140  i=l,  nt2 
x(i)=0.0 
y(i)=0.0 

z(i)=^.5+(i-nt2)*1.0 
x(i+nt2)=0.0 
y(i+nt2)=0.0 
z(i+nt2)=0.54(i-l  )*  1 .0 
140  continue 

else 

do  142  i=l,  nt2 

x0)=0.0 

y(i)=0.0 

2(i)=-l.Of(i-nt2)*1.0 

x(i+nt2+l)=0.0 
y(i+nt2+l)=0.0 
z(i+nt2+l)=1.0+<i-l)*L0 
142  continue 

x(nt2+l)=0.0 

y(nt2+l)=0.0 

z(nt2+l)=0.0 

end  if 
xl=nt 
return 
end 
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c 


c — s  u  b:  d)deq> 

c -  cholesky  decomposition  :a=(IXdXu) 

parameter  (nuip=l(X),  mn3=mnp*3) 
complex  a(mn3,mn3),  aij,  ajj 
do40j=l4i 

jH'l 

if(jl.lt.2)goto25 

do20i=2jl 

aij=cmplx(0.0,0.0) 

dol0k=141 

10  aij=aij+a(k4)*a(kj)*a(kjt) 

a(ij)=(a(ij)-aij)/a(i,i) 

20  aO‘4)=a0d) 

25  ajj=  cmplx(0.0,0.0) 

if  (j.eq.l)go  to  35 
do  30  kl=l jl 

30  ajj=ajj+  a(kl  j)*a(kl  j)*a(kl  Jcl ) 

aOJ)=aOd>aij 

35  if(a(j  j).eq.O.)  go  to  50 

40  continue 
return 
50  print  55 

55  fonnat(/  ’a(j  j)  is  zero  program  stop.') 
return 
end 
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c— sub:  coord  transforms  (x',y',z')  into  measuring  coord.(x,y^) 
couples  with  size  parameter  alpha*2 
subroutine  coord(np,aIpha,cchi,schi,cpsi,spsi) 
parameter  (mnp=100) 
dimension  c(33),p(3) 
common/1  l/xo(mnp),yo(mnp)^mnp) 
common/12/x(mnp),y(mnp),z(mnp) 
common/transl/c 
c(l,l>=cpsi*cchi 
c(l^)=-spsi 
c(l,3>=cpsi*schi 
c(2,l)=spsi*cchi 
c(2^)=cpsi 
c(2,3)=spsi*schi 
c(3,l)=-schi 

coa)=o. 

c(33)=cchi 
do  10  i=:l4>p 
p(l)=xo(i) 
p(2)=yoCi) 
p(3)=zo(i) 
x(i)=0. 
y(i)=0. 
zCi)=0. 
do  10  j=13 

x(i)=x(i>+2.*alpha*c(IJ)*P0) 

y(i)=y(i)+2.*alpha*c(2j)*p(j) 

10  z(i)=za}+2.*alpha*c(3J)*p(j) 
return 
end 
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0- 

c — subrcbsolx 

c - OX<JXuXx)=(b)with  (d)  cannot  be  zero 

c - let(lXz)=(b),and(dXuXx)=(z) 

c — ^verified  November  15, 1993 — 
subroutine  cbsoIx(a,b,n) 
parameter(mnp=100,  nin3=mnp*3) 
complex  a(mn34nn3),b(n) 
do70i=2,n 
il=i-l 

do60k=141 

60  b(i)=b{i)-a(k4)*b(k) 

70  continue 
c  divided  by  diagonal 
c  (dXyMz),  where  (y)  to  be  found 
do  72  i=l,n 
72  bCi)=b(i)/a(i4) 
c  backward  substitution 

c  (u)(x),  where  (x)  is  final  solution 

do80i  =n-l,l,-l 
do75  j=i+l,n 
75  b(i)=b(i)-a(i  j)*bO) 

80  continue 
return 
end 
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VIII.  USER'S  GUIDE  FOR  AGGL  WITH  SOURCE  CODE 


1.  User  Interface:  AGGL_Shell 

To  use  AGGL,  simply  call  AGGL_SheII  from  your  main  program  in  the  following  format: 

CALL  AGGL^Shell(W,NX,NPS  AJ^PHA,RFIN^ 

where 

ALPHA:  Primaiy  particle  size  parameter 

NP:  Number  of  primary  particles  in  the  agglonwrate 

NX:  Number  of  scattering  angles 

NPS :  Number  of  structures  used  to  mimic  randomly  branched  chains 

RFINDX:  Complex  refractive  index,  the  imaginaiy  part  must  be  positive 

Other  input/output  dataare  transferred  via  COMMON  blocks,  see  the  following  list. 

Block  Variables  Explanation 

name  carrying 

LI  X,Y,Z:  Coordinate  components  of  the  particle  centers  of  the  agglomerate 

P2  CHAIN  Type  of  chains,  ’R’  for  randomly  branched  chain,  'S’  for  straight  chain, 

’C  for  cluster,  'U'  for  user  specified  geometry 

LOU  LUl  Switch  for  Qsca  computation,  0  for  subtraction,  1  for  direct 

computation 

LU3  Switch  for  orientation,  0  for  random  orientaticm 

1  for  a  specific  orientation 

LU4  Switch  for  derivative  computation,  0  for  no  derivative  computation, 

1  to  calculate  derivatives 

Qextm  QE,QA,QS  Extinction,  absorption,  and  total  scattering  efficiency  factors 
Louv  RPP  Etissymmetry  ,  depolarization  ratios  and  extinction  normalized 

differential  scattering  efficiencies  under  vertical  inddence 
Louh  Same  as  Louv  but  under  horizontal  inddence 

J2sm  QPP  Differential  scattering  efficiencies  under  both  vertical  and  horizcmtal 

incidences.  P  can  be  V  or  H  standing  for  vertical  or  horizmital 
ANGK  Scattering  angles  where  diffierential  scattering  is  to  be  con^iuted 
LON  DPPN  Derivatives  with  respect  to  the  real  part  of  the  refractive  index 

LOK  DPPK  Derivatives  with  respect  to  the  imaginaiy  part  of  the  refractive  index 

Loux  DPPX  Derivatives  with  respect  to  the  size  parameter  x 
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Orient  Dchi,Dpsi  Orientation  of  the  agglomerate  if  not  iand(m 

2.  Output  List 

QEXT :  Extinction  efficiency  factor  for  vertical  incidence 

QABS:  Absorption  efficiency  factor  for  vertical  ind(knce 

QSCA:  Total  scattering  efficiency  factor  for  vertical  inddence 

QEXTH:  Same  as  QEXT,  but  for  horizontal  inddence 

QABSH:  Same  as  QABS,  but  for  horizontal  inddence 

QSCAH:  Same  as  QSCA,  but  for  horizontal  inddence 

CPP():  Differential  scattering  efficiencies.  P  can  be  either  V  or  H 

standing  for  vertical  or  horizontal  polarization  states.  The  first  P 

is  related  to  the  observation  or  measurement,  and  the  second  P  the  incidence  light. 

RPP():  EHssymmetry  or  depolarization  ratios(optional) 

KWQ:  Extinction  normalized  differential  scattering(optional) 

DQEXTq:  Derivative  of  QEXT  with  respect  to  q.  Q 

can  be  n(real  part  of  the  refractive  index),  k  (imaginary  part  of  the  refractive  index), 
or  X  (size  parameter) 

DCPPq:  Derivative  of  CPP  with  respect  to  q.  (Optional) 
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3  Sample  Master  Program  and  Results 

The  following  is  a  sample  master  program  which  makes  use  of  AGGL  to  calculate  the  extinction,  absorption, 
total  scattering,  and  differential  scattering  efficiencies  under  both  vertical  incidence  and  horizontal  incidence  for  a 
straight  chain  of  10  particles  with  a  specific  orientation  ^=45  deg  and  xp=60  deg).  Other  parameters  are  x=0.2  and 
m=L7+0.7i.  Computational  results  are  given  after  the  program. 

Note  the  symbols  ”$"  and  are  exclusively  used  to  mark  continuing  lines.  Other  symbols  used  for  the  same 
purpose  are  and 
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C —  THIS  IS  A  SAMPLE  USER'  MASTER  PROGRAM  **♦*♦***♦ 

C23456789012345678901234567890123456789012345678901234567890123456789 
C  PROGRAM  MASTER.FOR 

PARAMETER  (NS=181,  NC=1,  MNP^I00^=37) 

C***THE  FOLLOWING  LINES  CAN  NOT  BE  CHANGED.  FOLLOW  STRICTLY!*** 

C***  WHILE  WRITING  YOUR  OWN  MASTER  PROGRAM  ****** 

COMPLEX  RFINDX,EPSLON 
CHARACTER  EFF*I,POL*l,CHAIN*l 
COMMON/P2/POL,EFF,CHAIN 
COMMON/LOU/LU,LUl  T-U2T.U3,LU4 

COMMON/J2SM/SCAHHK(181),SCAWK(181),  SCAHVK(I81),  SCAVHK(181), 

$  ANG(181) 

COMMON/QEXTM/QEXT,QABS,QSCA,QEXTH,QABSH,QSCAH 

COMMON/ORIENT/DCHI,DPSI 

REAL  KVVE(NS)4CHHE(NS).KHVE(NS),KVHE(NS) 

COMMON/LOUV/RW(NS)T)RVVN(NS),DRWK(NS), 

$  DKHVEN(NS)J)KHVEK(NS), 

$  KVVE,DKWEN(NS)4)KWEK(NS),RHV(NS)T)RHVN(NS), 

$  DRHVK(NS),KHVE,KVHE 
COMMON/LOUH/RHH(NS),DRHHN(NS),DRHHK(NS), 

$  DKVHEN(NS),DKVHEK(NS), 

$  KHHE,DKHHEN(NS)JDKHHEK(NS)TIVH(NS)T)RVHN(NS)4)RVHK(NS) 
COMMON/LOUX/QETXT)RWX(NSXDRHVX(NS),DRVHX(NS),DRHHX(NS), 

$  QEHTX,DKWEX(NS)T)KHVEX(NS),DKVHEX(NS)4)KHHEX(NS) 

C***ALL  OF  THE  NAMES  APPEARED  ABOVE  ARE  RESERVED.  *************** 
C***  DO  NOT  ATTEMPT  TO  CHANGE  OR  RENAME  *********** 

CHARACTER*40  NAMDAT  J4AMOUTJ>IAMPLT 
NAM0UT='[MEL0U.MIE]CHAIS6.0UT’ 

NAMDAT='fMELOU.MIE]CHAIS6.DAr 

NAMPLT='[MELOU.MIE]CHAS.DAT' 

OPEN(UNrT=12J=ILE=NAMOUT  ,STATUS=’UNKNOWN') 
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0PEN(UNIT=133LE=NAMPLT,STATUS='UNKN0WN') 

OPEN(UMT=143LE=NANroAT,STATUS='UNKNOWN') 


c=opnc»sis  . . . :■  -  . 

c— SELF-TERM  OPTIONS:  LU 
C  LU=0  FOR  JONES; 

C—  1  FOR  THIRD  ORDER  APPROXIMATION  (GOEDBCKE«&0'BRIEN) 

C—  2  FOR  EXACT  SOLUTION 

C- - 3  OTHERS 

LU=1 


C— TOTAL  SCATTERING  OPTIONS: 

C  LU1=0  FOR  SUBSTRACTION  FROM  EXTINCTION  AND  ABSORPHON 

C - 1  FOR  DIRECT  INTEGRATION 


LU1=0 


C— T  MATRIX  OPTIONS: 

C  LU2^  FOR  JONES  (INCORRECTVDISABLED 

C - 1  FOR  CORRECTED  TENSOR  MATRIX 


LU2=1 

C--ORIENTATION  OPTIONS: 

C  LU3=0  FOR  RANDOM  ORIENTATION; 

C -  1  FOR  A  SPECMC  ORIENTATION 

LU3=1 

IF(LU3.EQ.1)THEN 
DCHI=45.0 
DPSI=60.0 
END  IF 

C— SCATTERING  DERIVATIVES  OPTIONS: 
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C  LU4=1  FOR  DERIVATIVE  CALCULATIONS 

C - 0  TO  SKIP  DERIVATIVES 

LU4=0 

C— OTHER  OPTIONS/  POL  OPTION  DISABLED/ 

C— CHAIN:  SELECT  AMONG  ’S’,  'C,  'R',  OR,  "U' 

EFF='P 

POI^H’ 

CHAIN='S’ 

C— DATA  SPEanCATION; 

IF(CHAIN.EQ.’R')THEN 
NPS=5{)0 
END  IF 
NP=  10 
DO  100  K=1,NX 

100  ANG(K)=(K-1)*180.0/(NX-1) 

WL=.488 
REF  =1.7 
RIM  =.7 
PI=3. 1415926 
REFMED=1.0 

RFINDX=CMPLX(REF,  RIM)/REFMED 

DO303IL=l,NC 

ALPHA=0.2+0.05*(IL-1) 

DIA=ALPHA*WL/PI 

CALL  AGGL_SHELL(NPJDCNPS,ALPHA  JIHNDX) 

C!!!— DATA  OUTPUT  BLOCK - 

C!!!— THIS  BLOCK  IS  NOT  CRUTIAL  AND  OMITTED  HERE  TO  SAVE  SPACE 


303  CONTINUE 
STOP 
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END 


C***This  is  the  output  of  the  sample  master  program.  Make  sure  run  the  preceding  sample  C***program  and  obtain 
the  same  result  before  using  AGGL_Shell  for  your  applications 


DATAHLES: 


(MEL0U.inie]CHAIs6.0UT 


[MELOU.Mie]CHAIs6.DAT 

PRIMARY  PARTICLE  DIAMETER  =  0.03107  UM 

REFRACTIVE  INDEX:  M  =  1 .700  -1 0.700 

NUMBER  OF  PRIMARY  PARTICLES  =  10 

PRIMARY  PARTICLE  SIZE  PARAMETER  =  0.200 

AGGLOMERATE  SIZE  PARAMETER  =  2.000 

STRAIGHT  CHAINS 

FULLY  COUPLED  SOLUTION 

NUMBER  OF  ITERATIONS  1 

SPECIFIED  ORIENTATION 

CHI=  45.00000  DEGSAND  PSI=  60.00000 
SCATTERING  COEFHaENT(VERT.)  QSCAT=  6.5750182E-03 
EXTINCTION  COEFnaENT(VERT.)QEXT=  0.2265852 
ABSORPTION  COSnCIENT(VERT.)QABS=  0.2200102 
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SCATTERING  COEFnaENT(Hort.)  QSCATh; 
EXTINCTION  COEFnCIENT(Hort)  QEXTh= 
ABSORPTION  COEFFICIENT(HorL)  QABSh= 


=  8.2514584E-03 

0.2679992 

Q259rJ411 
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Scattering 

Angle  (deg)  Qyy 

0.0  0.1468E-02 

5.0  0.1463E-02 

10.0  0.1446E-02 

15.0  0.1420E-02 

20.0  0.1390E-02 

25.0  0.1358E-02 

30.0  0.1327E-02 

35.0  0.1300E-02 

40.0  0.1277E-02 

45.0  0.1261E-02 

50.0  0.1252E-02 

55.0  0.1249E-02 

60.0  0.1253E-02 

65.0  0.1263E-02 

70.0  0.1278E-02 

75.0  0.1295E-02 

80.0  0.1313E-02 

85.0  0.1329E-02 

90.0  0.1340E-02 

95.0  0.1343E-02 

100.0  0.1336E-02 

105.0  0.1316E-02 

110.0  0.1281E-02 

115.0  0.1230E-02 

120.0  0.1165E-02 

125.0  0.1086E-02 

130.0  0.9966E-03 

135.0  0.8991E-03 

140.0  0.7975E-03 

145.0  0.6954E-03 

150.0  0.5962E-03 

155.0  0.5029E-03 


Differential  Scattering  Efficiencies 

Qhh  Qhv  Qv 


0.1709E-02 

0.1663E-02 

0.1581E-02 

0.1471E-02 

0.1340E-02 

0.1197E-02 

0.1048E-02 

0.8992E-03 

0.7545E-03 

0.6171E-03 

0.4894E-03 

0.3730E-03 

0.2694E-03 

0.1802E-03 

0.1072E-03 

0.5243E-04 

0.1822E-04 

0.6583E-05 

0.18%E-04 

0.5562E-04 

0.1151E-03 

0.1939E-03 

0.2864E-03 

0.3852E-03 

0.4816E-03 

0.5672E-03 

0.6343E-03 

0.6775E-03 

0.6940E-03 

0.6840E-03 

0.6504E-03 

0.5980E-03 


0.1128E-04 

0.9565E-05 

0.7863E-05 

0.6242E-05 

0.4755E-05 

0.3440E-05 

0.2323E-05 

0.1418E-05 

0.7332E-06 

0.2715E-06 

0.3435E-07 

0.2302E-07 

0.2397E-06 

0.6871E-06 

0.1368E-05 

0.2281E-05 

0.3423E-05 

0.4778E-05 

0.6322E-05 

0.8013E-05 

0.9797E-05 

0.1160E-04 

0.1335E-04 

0.1495E-04 

0.1632E-04 

0.1739E-04 

0.1808E-04 

0.1838E-04 

0.1827E-04 

0.1778E-04 

0.1695E-04 

0.1585E-04 


0.1128E-04 

O.llllE-04 

0.1090E-04 

0.1066E-04 

0.1041E4)4 

0.1015E-04 

0.9905E-05 

0.9675E-05 

0.9472E-05 

0.9303E-05 

0.9173E-05 

0.9087E-05 

0.9047E-05 

0.9056E-05 

0.9114E-05 

0.9219E-05 

0.9369E-05 

0.9559E-05 

0.9784E-05 

0.1004E-04 

0.1030E-04 

0.1058E-04 

0.1084E-04 

0.1109E-04 

0.1130E-04 

0.1146E-04 

0.1156E-04 

0.1159E-04 

0.1154E-04 

0.1141E-04 

0.1120E-04 

0.1092E-04 
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160.0 

0.4175E-03 

0.5326E-03 

0.1454E-04 

0.1057E-04 

165.0 

0.3413E-03 

0.4603E-03 

0.1310E-04 

0.1015E-04 

170.0 

0.2750E-03 

0.3868E-03 

0.1161E-04 

0.%99E-05 

175.0 

0.2186E-03 

0.3163E-03 

0.1013E04 

0.9212E-05 

180.0 

0.1715E-03 

0.2521E-03 

0.8709E-05 

0.8709E-05 
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